home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / zgesvd.f < prev    next >
Text File  |  1996-07-19  |  142KB  |  3,611 lines
  1.       SUBROUTINE ZGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,
  2.      $                   WORK, LWORK, RWORK, INFO )
  3. *
  4. *  -- LAPACK driver routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       CHARACTER          JOBU, JOBVT
  11.       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
  12. *     ..
  13. *     .. Array Arguments ..
  14.       DOUBLE PRECISION   RWORK( * ), S( * )
  15.       COMPLEX*16         A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
  16.      $                   WORK( * )
  17. *     ..
  18. *
  19. *  Purpose
  20. *  =======
  21. *
  22. *  ZGESVD computes the singular value decomposition (SVD) of a complex
  23. *  M-by-N matrix A, optionally computing the left and/or right singular
  24. *  vectors. The SVD is written
  25. *
  26. *       A = U * SIGMA * conjugate-transpose(V)
  27. *
  28. *  where SIGMA is an M-by-N matrix which is zero except for its
  29. *  min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
  30. *  V is an N-by-N unitary matrix.  The diagonal elements of SIGMA
  31. *  are the singular values of A; they are real and non-negative, and
  32. *  are returned in descending order.  The first min(m,n) columns of
  33. *  U and V are the left and right singular vectors of A.
  34. *
  35. *  Note that the routine returns V**H, not V.
  36. *
  37. *  Arguments
  38. *  =========
  39. *
  40. *  JOBU    (input) CHARACTER*1
  41. *          Specifies options for computing all or part of the matrix U:
  42. *          = 'A':  all M columns of U are returned in array U:
  43. *          = 'S':  the first min(m,n) columns of U (the left singular
  44. *                  vectors) are returned in the array U;
  45. *          = 'O':  the first min(m,n) columns of U (the left singular
  46. *                  vectors) are overwritten on the array A;
  47. *          = 'N':  no columns of U (no left singular vectors) are
  48. *                  computed.
  49. *
  50. *  JOBVT   (input) CHARACTER*1
  51. *          Specifies options for computing all or part of the matrix
  52. *          V**H:
  53. *          = 'A':  all N rows of V**H are returned in the array VT;
  54. *          = 'S':  the first min(m,n) rows of V**H (the right singular
  55. *                  vectors) are returned in the array VT;
  56. *          = 'O':  the first min(m,n) rows of V**H (the right singular
  57. *                  vectors) are overwritten on the array A;
  58. *          = 'N':  no rows of V**H (no right singular vectors) are
  59. *                  computed.
  60. *
  61. *          JOBVT and JOBU cannot both be 'O'.
  62. *
  63. *  M       (input) INTEGER
  64. *          The number of rows of the input matrix A.  M >= 0.
  65. *
  66. *  N       (input) INTEGER
  67. *          The number of columns of the input matrix A.  N >= 0.
  68. *
  69. *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
  70. *          On entry, the M-by-N matrix A.
  71. *          On exit,
  72. *          if JOBU = 'O',  A is overwritten with the first min(m,n)
  73. *                          columns of U (the left singular vectors,
  74. *                          stored columnwise);
  75. *          if JOBVT = 'O', A is overwritten with the first min(m,n)
  76. *                          rows of V**H (the right singular vectors,
  77. *                          stored rowwise);
  78. *          if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
  79. *                          are destroyed.
  80. *
  81. *  LDA     (input) INTEGER
  82. *          The leading dimension of the array A.  LDA >= max(1,M).
  83. *
  84. *  S       (output) DOUBLE PRECISION array, dimension (min(M,N))
  85. *          The singular values of A, sorted so that S(i) >= S(i+1).
  86. *
  87. *  U       (output) COMPLEX*16 array, dimension (LDU,UCOL)
  88. *          (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
  89. *          If JOBU = 'A', U contains the M-by-M unitary matrix U;
  90. *          if JOBU = 'S', U contains the first min(m,n) columns of U
  91. *          (the left singular vectors, stored columnwise);
  92. *          if JOBU = 'N' or 'O', U is not referenced.
  93. *
  94. *  LDU     (input) INTEGER
  95. *          The leading dimension of the array U.  LDU >= 1; if
  96. *          JOBU = 'S' or 'A', LDU >= M.
  97. *
  98. *  VT      (output) COMPLEX*16 array, dimension (LDVT,N)
  99. *          If JOBVT = 'A', VT contains the N-by-N unitary matrix
  100. *          V**H;
  101. *          if JOBVT = 'S', VT contains the first min(m,n) rows of
  102. *          V**H (the right singular vectors, stored rowwise);
  103. *          if JOBVT = 'N' or 'O', VT is not referenced.
  104. *
  105. *  LDVT    (input) INTEGER
  106. *          The leading dimension of the array VT.  LDVT >= 1; if
  107. *          JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
  108. *
  109. *  WORK    (workspace/output) COMPLEX*16 array, dimension (LWORK)
  110. *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  111. *
  112. *  LWORK   (input) INTEGER
  113. *          The dimension of the array WORK. LWORK >= 1.
  114. *          LWORK >=  2*MIN(M,N)+MAX(M,N).
  115. *          For good performance, LWORK should generally be larger.
  116. *
  117. *  RWORK   (workspace) DOUBLE PRECISION array, dimension
  118. *                                  (max(3*min(M,N),5*min(M,N)-4))
  119. *          On exit, if INFO > 0, RWORK(1:MIN(M,N)-1) contains the
  120. *          unconverged superdiagonal elements of an upper bidiagonal
  121. *          matrix B whose diagonal is in S (not necessarily sorted).
  122. *          B satisfies A = U * B * VT, so it has the same singular
  123. *          values as A, and singular vectors related by U and VT.
  124. *
  125. *  INFO    (output) INTEGER
  126. *          = 0:  successful exit.
  127. *          < 0:  if INFO = -i, the i-th argument had an illegal value.
  128. *          > 0:  if ZBDSQR did not converge, INFO specifies how many
  129. *                superdiagonals of an intermediate bidiagonal form B
  130. *                did not converge to zero. See the description of RWORK
  131. *                above for details.
  132. *
  133. *  =====================================================================
  134. *
  135. *     .. Parameters ..
  136.       COMPLEX*16         CZERO, CONE
  137.       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
  138.      $                   CONE = ( 1.0D0, 0.0D0 ) )
  139.       DOUBLE PRECISION   ZERO, ONE
  140.       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  141. *     ..
  142. *     .. Local Scalars ..
  143.       LOGICAL            WNTUA, WNTUAS, WNTUN, WNTUO, WNTUS, WNTVA,
  144.      $                   WNTVAS, WNTVN, WNTVO, WNTVS
  145.       INTEGER            BLK, CHUNK, I, IE, IERR, IR, IRWORK, ISCL,
  146.      $                   ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU,
  147.      $                   MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU,
  148.      $                   NRVT, WRKBL
  149.       DOUBLE PRECISION   ANRM, BIGNUM, EPS, SMLNUM
  150. *     ..
  151. *     .. Local Arrays ..
  152.       DOUBLE PRECISION   DUM( 1 )
  153.       COMPLEX*16         CDUM( 1 )
  154. *     ..
  155. *     .. External Subroutines ..
  156.       EXTERNAL           DLASCL, XERBLA, ZBDSQR, ZGEBRD, ZGELQF, ZGEMM,
  157.      $                   ZGEQRF, ZLACPY, ZLASCL, ZLASET, ZUNGBR, ZUNGLQ,
  158.      $                   ZUNGQR, ZUNMBR
  159. *     ..
  160. *     .. External Functions ..
  161.       LOGICAL            LSAME
  162.       INTEGER            ILAENV
  163.       DOUBLE PRECISION   DLAMCH, ZLANGE
  164.       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE
  165. *     ..
  166. *     .. Intrinsic Functions ..
  167.       INTRINSIC          MAX, MIN, SQRT
  168. *     ..
  169. *     .. Executable Statements ..
  170. *
  171. *     Test the input arguments
  172. *
  173.       INFO = 0
  174.       MINMN = MIN( M, N )
  175.       MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 )
  176.       WNTUA = LSAME( JOBU, 'A' )
  177.       WNTUS = LSAME( JOBU, 'S' )
  178.       WNTUAS = WNTUA .OR. WNTUS
  179.       WNTUO = LSAME( JOBU, 'O' )
  180.       WNTUN = LSAME( JOBU, 'N' )
  181.       WNTVA = LSAME( JOBVT, 'A' )
  182.       WNTVS = LSAME( JOBVT, 'S' )
  183.       WNTVAS = WNTVA .OR. WNTVS
  184.       WNTVO = LSAME( JOBVT, 'O' )
  185.       WNTVN = LSAME( JOBVT, 'N' )
  186.       MINWRK = 1
  187. *
  188.       IF( .NOT.( WNTUA .OR. WNTUS .OR. WNTUO .OR. WNTUN ) ) THEN
  189.          INFO = -1
  190.       ELSE IF( .NOT.( WNTVA .OR. WNTVS .OR. WNTVO .OR. WNTVN ) .OR.
  191.      $         ( WNTVO .AND. WNTUO ) ) THEN
  192.          INFO = -2
  193.       ELSE IF( M.LT.0 ) THEN
  194.          INFO = -3
  195.       ELSE IF( N.LT.0 ) THEN
  196.          INFO = -4
  197.       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  198.          INFO = -6
  199.       ELSE IF( LDU.LT.1 .OR. ( WNTUAS .AND. LDU.LT.M ) ) THEN
  200.          INFO = -9
  201.       ELSE IF( LDVT.LT.1 .OR. ( WNTVA .AND. LDVT.LT.N ) .OR.
  202.      $         ( WNTVS .AND. LDVT.LT.MINMN ) ) THEN
  203.          INFO = -11
  204.       END IF
  205. *
  206. *     Compute workspace
  207. *      (Note: Comments in the code beginning "Workspace:" describe the
  208. *       minimal amount of workspace needed at that point in the code,
  209. *       as well as the preferred amount for good performance.
  210. *       CWorkspace refers to complex workspace, and RWorkspace to
  211. *       real workspace. NB refers to the optimal block size for the
  212. *       immediately following subroutine, as returned by ILAENV.)
  213. *
  214.       IF( INFO.EQ.0 .AND. LWORK.GE.1 .AND. M.GT.0 .AND. N.GT.0 ) THEN
  215.          IF( M.GE.N ) THEN
  216. *
  217. *           Space needed for ZBDSQR is BDSPAC = MAX( 3*N, 5*N-4 )
  218. *
  219.             IF( M.GE.MNTHR ) THEN
  220.                IF( WNTUN ) THEN
  221. *
  222. *                 Path 1 (M much larger than N, JOBU='N')
  223. *
  224.                   MAXWRK = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1,
  225.      $                     -1 )
  226.                   MAXWRK = MAX( MAXWRK, 2*N+2*N*
  227.      $                     ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
  228.                   IF( WNTVO .OR. WNTVAS )
  229.      $               MAXWRK = MAX( MAXWRK, 2*N+( N-1 )*
  230.      $                        ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )
  231.                   MINWRK = 3*N
  232.                   MAXWRK = MAX( MINWRK, MAXWRK )
  233.                ELSE IF( WNTUO .AND. WNTVN ) THEN
  234. *
  235. *                 Path 2 (M much larger than N, JOBU='O', JOBVT='N')
  236. *
  237.                   WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
  238.                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M,
  239.      $                    N, N, -1 ) )
  240.                   WRKBL = MAX( WRKBL, 2*N+2*N*
  241.      $                    ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
  242.                   WRKBL = MAX( WRKBL, 2*N+N*
  243.      $                    ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) )
  244.                   MAXWRK = MAX( N*N+WRKBL, N*N+M*N )
  245.                   MINWRK = 2*N + M
  246.                   MAXWRK = MAX( MINWRK, MAXWRK )
  247.                ELSE IF( WNTUO .AND. WNTVAS ) THEN
  248. *
  249. *                 Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
  250. *                 'A')
  251. *
  252.                   WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
  253.                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M,
  254.      $                    N, N, -1 ) )
  255.                   WRKBL = MAX( WRKBL, 2*N+2*N*
  256.      $                    ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
  257.                   WRKBL = MAX( WRKBL, 2*N+N*
  258.      $                    ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) )
  259.                   WRKBL = MAX( WRKBL, 2*N+( N-1 )*
  260.      $                    ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )
  261.                   MAXWRK = MAX( N*N+WRKBL, N*N+M*N )
  262.                   MINWRK = 2*N + M
  263.                   MAXWRK = MAX( MINWRK, MAXWRK )
  264.                ELSE IF( WNTUS .AND. WNTVN ) THEN
  265. *
  266. *                 Path 4 (M much larger than N, JOBU='S', JOBVT='N')
  267. *
  268.                   WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
  269.                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M,
  270.      $                    N, N, -1 ) )
  271.                   WRKBL = MAX( WRKBL, 2*N+2*N*
  272.      $                    ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
  273.                   WRKBL = MAX( WRKBL, 2*N+N*
  274.      $                    ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) )
  275.                   MAXWRK = N*N + WRKBL
  276.                   MINWRK = 2*N + M
  277.                   MAXWRK = MAX( MINWRK, MAXWRK )
  278.                ELSE IF( WNTUS .AND. WNTVO ) THEN
  279. *
  280. *                 Path 5 (M much larger than N, JOBU='S', JOBVT='O')
  281. *
  282.                   WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
  283.                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M,
  284.      $                    N, N, -1 ) )
  285.                   WRKBL = MAX( WRKBL, 2*N+2*N*
  286.      $                    ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
  287.                   WRKBL = MAX( WRKBL, 2*N+N*
  288.      $                    ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) )
  289.                   WRKBL = MAX( WRKBL, 2*N+( N-1 )*
  290.      $                    ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )
  291.                   MAXWRK = 2*N*N + WRKBL
  292.                   MINWRK = 2*N + M
  293.                   MAXWRK = MAX( MINWRK, MAXWRK )
  294.                ELSE IF( WNTUS .AND. WNTVAS ) THEN
  295. *
  296. *                 Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
  297. *                 'A')
  298. *
  299.                   WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
  300.                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M,
  301.      $                    N, N, -1 ) )
  302.                   WRKBL = MAX( WRKBL, 2*N+2*N*
  303.      $                    ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
  304.                   WRKBL = MAX( WRKBL, 2*N+N*
  305.      $                    ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) )
  306.                   WRKBL = MAX( WRKBL, 2*N+( N-1 )*
  307.      $                    ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )
  308.                   MAXWRK = N*N + WRKBL
  309.                   MINWRK = 2*N + M
  310.                   MAXWRK = MAX( MINWRK, MAXWRK )
  311.                ELSE IF( WNTUA .AND. WNTVN ) THEN
  312. *
  313. *                 Path 7 (M much larger than N, JOBU='A', JOBVT='N')
  314. *
  315.                   WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
  316.                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'ZUNGQR', ' ', M,
  317.      $                    M, N, -1 ) )
  318.                   WRKBL = MAX( WRKBL, 2*N+2*N*
  319.      $                    ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
  320.                   WRKBL = MAX( WRKBL, 2*N+N*
  321.      $                    ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) )
  322.                   MAXWRK = N*N + WRKBL
  323.                   MINWRK = 2*N + M
  324.                   MAXWRK = MAX( MINWRK, MAXWRK )
  325.                ELSE IF( WNTUA .AND. WNTVO ) THEN
  326. *
  327. *                 Path 8 (M much larger than N, JOBU='A', JOBVT='O')
  328. *
  329.                   WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
  330.                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'ZUNGQR', ' ', M,
  331.      $                    M, N, -1 ) )
  332.                   WRKBL = MAX( WRKBL, 2*N+2*N*
  333.      $                    ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
  334.                   WRKBL = MAX( WRKBL, 2*N+N*
  335.      $                    ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) )
  336.                   WRKBL = MAX( WRKBL, 2*N+( N-1 )*
  337.      $                    ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )
  338.                   MAXWRK = 2*N*N + WRKBL
  339.                   MINWRK = 2*N + M
  340.                   MAXWRK = MAX( MINWRK, MAXWRK )
  341.                ELSE IF( WNTUA .AND. WNTVAS ) THEN
  342. *
  343. *                 Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
  344. *                 'A')
  345. *
  346.                   WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
  347.                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'ZUNGQR', ' ', M,
  348.      $                    M, N, -1 ) )
  349.                   WRKBL = MAX( WRKBL, 2*N+2*N*
  350.      $                    ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
  351.                   WRKBL = MAX( WRKBL, 2*N+N*
  352.      $                    ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) )
  353.                   WRKBL = MAX( WRKBL, 2*N+( N-1 )*
  354.      $                    ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )
  355.                   MAXWRK = N*N + WRKBL
  356.                   MINWRK = 2*N + M
  357.                   MAXWRK = MAX( MINWRK, MAXWRK )
  358.                END IF
  359.             ELSE
  360. *
  361. *              Path 10 (M at least N, but not much larger)
  362. *
  363.                MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N,
  364.      $                  -1, -1 )
  365.                IF( WNTUS .OR. WNTUO )
  366.      $            MAXWRK = MAX( MAXWRK, 2*N+N*
  367.      $                     ILAENV( 1, 'ZUNGBR', 'Q', M, N, N, -1 ) )
  368.                IF( WNTUA )
  369.      $            MAXWRK = MAX( MAXWRK, 2*N+M*
  370.      $                     ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) )
  371.                IF( .NOT.WNTVN )
  372.      $            MAXWRK = MAX( MAXWRK, 2*N+( N-1 )*
  373.      $                     ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )
  374.                MINWRK = 2*N + M
  375.                MAXWRK = MAX( MINWRK, MAXWRK )
  376.             END IF
  377.          ELSE
  378. *
  379. *           Space needed for ZBDSQR is BDSPAC = MAX( 3*M, 5*M-4 )
  380. *
  381.             IF( N.GE.MNTHR ) THEN
  382.                IF( WNTVN ) THEN
  383. *
  384. *                 Path 1t(N much larger than M, JOBVT='N')
  385. *
  386.                   MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1,
  387.      $                     -1 )
  388.                   MAXWRK = MAX( MAXWRK, 2*M+2*M*
  389.      $                     ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
  390.                   IF( WNTUO .OR. WNTUAS )
  391.      $               MAXWRK = MAX( MAXWRK, 2*M+M*
  392.      $                        ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) )
  393.                   MINWRK = 3*M
  394.                   MAXWRK = MAX( MINWRK, MAXWRK )
  395.                ELSE IF( WNTVO .AND. WNTUN ) THEN
  396. *
  397. *                 Path 2t(N much larger than M, JOBU='N', JOBVT='O')
  398. *
  399.                   WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
  400.                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M,
  401.      $                    N, M, -1 ) )
  402.                   WRKBL = MAX( WRKBL, 2*M+2*M*
  403.      $                    ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
  404.                   WRKBL = MAX( WRKBL, 2*M+( M-1 )*
  405.      $                    ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) )
  406.                   MAXWRK = MAX( M*M+WRKBL, M*M+M*N )
  407.                   MINWRK = 2*M + N
  408.                   MAXWRK = MAX( MINWRK, MAXWRK )
  409.                ELSE IF( WNTVO .AND. WNTUAS ) THEN
  410. *
  411. *                 Path 3t(N much larger than M, JOBU='S' or 'A',
  412. *                 JOBVT='O')
  413. *
  414.                   WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
  415.                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M,
  416.      $                    N, M, -1 ) )
  417.                   WRKBL = MAX( WRKBL, 2*M+2*M*
  418.      $                    ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
  419.                   WRKBL = MAX( WRKBL, 2*M+( M-1 )*
  420.      $                    ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) )
  421.                   WRKBL = MAX( WRKBL, 2*M+M*
  422.      $                    ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) )
  423.                   MAXWRK = MAX( M*M+WRKBL, M*M+M*N )
  424.                   MINWRK = 2*M + N
  425.                   MAXWRK = MAX( MINWRK, MAXWRK )
  426.                ELSE IF( WNTVS .AND. WNTUN ) THEN
  427. *
  428. *                 Path 4t(N much larger than M, JOBU='N', JOBVT='S')
  429. *
  430.                   WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
  431.                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M,
  432.      $                    N, M, -1 ) )
  433.                   WRKBL = MAX( WRKBL, 2*M+2*M*
  434.      $                    ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
  435.                   WRKBL = MAX( WRKBL, 2*M+( M-1 )*
  436.      $                    ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) )
  437.                   MAXWRK = M*M + WRKBL
  438.                   MINWRK = 2*M + N
  439.                   MAXWRK = MAX( MINWRK, MAXWRK )
  440.                ELSE IF( WNTVS .AND. WNTUO ) THEN
  441. *
  442. *                 Path 5t(N much larger than M, JOBU='O', JOBVT='S')
  443. *
  444.                   WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
  445.                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M,
  446.      $                    N, M, -1 ) )
  447.                   WRKBL = MAX( WRKBL, 2*M+2*M*
  448.      $                    ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
  449.                   WRKBL = MAX( WRKBL, 2*M+( M-1 )*
  450.      $                    ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) )
  451.                   WRKBL = MAX( WRKBL, 2*M+M*
  452.      $                    ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) )
  453.                   MAXWRK = 2*M*M + WRKBL
  454.                   MINWRK = 2*M + N
  455.                   MAXWRK = MAX( MINWRK, MAXWRK )
  456.                ELSE IF( WNTVS .AND. WNTUAS ) THEN
  457. *
  458. *                 Path 6t(N much larger than M, JOBU='S' or 'A',
  459. *                 JOBVT='S')
  460. *
  461.                   WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
  462.                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M,
  463.      $                    N, M, -1 ) )
  464.                   WRKBL = MAX( WRKBL, 2*M+2*M*
  465.      $                    ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
  466.                   WRKBL = MAX( WRKBL, 2*M+( M-1 )*
  467.      $                    ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) )
  468.                   WRKBL = MAX( WRKBL, 2*M+M*
  469.      $                    ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) )
  470.                   MAXWRK = M*M + WRKBL
  471.                   MINWRK = 2*M + N
  472.                   MAXWRK = MAX( MINWRK, MAXWRK )
  473.                ELSE IF( WNTVA .AND. WNTUN ) THEN
  474. *
  475. *                 Path 7t(N much larger than M, JOBU='N', JOBVT='A')
  476. *
  477.                   WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
  478.                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'ZUNGLQ', ' ', N,
  479.      $                    N, M, -1 ) )
  480.                   WRKBL = MAX( WRKBL, 2*M+2*M*
  481.      $                    ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
  482.                   WRKBL = MAX( WRKBL, 2*M+( M-1 )*
  483.      $                    ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) )
  484.                   MAXWRK = M*M + WRKBL
  485.                   MINWRK = 2*M + N
  486.                   MAXWRK = MAX( MINWRK, MAXWRK )
  487.                ELSE IF( WNTVA .AND. WNTUO ) THEN
  488. *
  489. *                 Path 8t(N much larger than M, JOBU='O', JOBVT='A')
  490. *
  491.                   WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
  492.                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'ZUNGLQ', ' ', N,
  493.      $                    N, M, -1 ) )
  494.                   WRKBL = MAX( WRKBL, 2*M+2*M*
  495.      $                    ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
  496.                   WRKBL = MAX( WRKBL, 2*M+( M-1 )*
  497.      $                    ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) )
  498.                   WRKBL = MAX( WRKBL, 2*M+M*
  499.      $                    ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) )
  500.                   MAXWRK = 2*M*M + WRKBL
  501.                   MINWRK = 2*M + N
  502.                   MAXWRK = MAX( MINWRK, MAXWRK )
  503.                ELSE IF( WNTVA .AND. WNTUAS ) THEN
  504. *
  505. *                 Path 9t(N much larger than M, JOBU='S' or 'A',
  506. *                 JOBVT='A')
  507. *
  508.                   WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
  509.                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'ZUNGLQ', ' ', N,
  510.      $                    N, M, -1 ) )
  511.                   WRKBL = MAX( WRKBL, 2*M+2*M*
  512.      $                    ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
  513.                   WRKBL = MAX( WRKBL, 2*M+( M-1 )*
  514.      $                    ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) )
  515.                   WRKBL = MAX( WRKBL, 2*M+M*
  516.      $                    ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) )
  517.                   MAXWRK = M*M + WRKBL
  518.                   MINWRK = 2*M + N
  519.                   MAXWRK = MAX( MINWRK, MAXWRK )
  520.                END IF
  521.             ELSE
  522. *
  523. *              Path 10t(N greater than M, but not much larger)
  524. *
  525.                MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N,
  526.      $                  -1, -1 )
  527.                IF( WNTVS .OR. WNTVO )
  528.      $            MAXWRK = MAX( MAXWRK, 2*M+M*
  529.      $                     ILAENV( 1, 'ZUNGBR', 'P', M, N, M, -1 ) )
  530.                IF( WNTVA )
  531.      $            MAXWRK = MAX( MAXWRK, 2*M+N*
  532.      $                     ILAENV( 1, 'ZUNGBR', 'P', N, N, M, -1 ) )
  533.                IF( .NOT.WNTUN )
  534.      $            MAXWRK = MAX( MAXWRK, 2*M+( M-1 )*
  535.      $                     ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) )
  536.                MINWRK = 2*M + N
  537.                MAXWRK = MAX( MINWRK, MAXWRK )
  538.             END IF
  539.          END IF
  540.          WORK( 1 ) = MAXWRK
  541.       END IF
  542. *
  543.       IF( LWORK.LT.MINWRK ) THEN
  544.          INFO = -13
  545.       END IF
  546.       IF( INFO.NE.0 ) THEN
  547.          CALL XERBLA( 'ZGESVD', -INFO )
  548.          RETURN
  549.       END IF
  550. *
  551. *     Quick return if possible
  552. *
  553.       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
  554.          IF( LWORK.GE.1 )
  555.      $      WORK( 1 ) = ONE
  556.          RETURN
  557.       END IF
  558. *
  559. *     Get machine constants
  560. *
  561.       EPS = DLAMCH( 'P' )
  562.       SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
  563.       BIGNUM = ONE / SMLNUM
  564. *
  565. *     Scale A if max element outside range [SMLNUM,BIGNUM]
  566. *
  567.       ANRM = ZLANGE( 'M', M, N, A, LDA, DUM )
  568.       ISCL = 0
  569.       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
  570.          ISCL = 1
  571.          CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
  572.       ELSE IF( ANRM.GT.BIGNUM ) THEN
  573.          ISCL = 1
  574.          CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
  575.       END IF
  576. *
  577.       IF( M.GE.N ) THEN
  578. *
  579. *        A has at least as many rows as columns. If A has sufficiently
  580. *        more rows than columns, first reduce using the QR
  581. *        decomposition (if sufficient workspace available)
  582. *
  583.          IF( M.GE.MNTHR ) THEN
  584. *
  585.             IF( WNTUN ) THEN
  586. *
  587. *              Path 1 (M much larger than N, JOBU='N')
  588. *              No left singular vectors to be computed
  589. *
  590.                ITAU = 1
  591.                IWORK = ITAU + N
  592. *
  593. *              Compute A=Q*R
  594. *              (CWorkspace: need 2*N, prefer N+N*NB)
  595. *              (RWorkspace: need 0)
  596. *
  597.                CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
  598.      $                      LWORK-IWORK+1, IERR )
  599. *
  600. *              Zero out below R
  601. *
  602.                CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
  603.      $                      LDA )
  604.                IE = 1
  605.                ITAUQ = 1
  606.                ITAUP = ITAUQ + N
  607.                IWORK = ITAUP + N
  608. *
  609. *              Bidiagonalize R in A
  610. *              (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
  611. *              (RWorkspace: need N)
  612. *
  613.                CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
  614.      $                      WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
  615.      $                      IERR )
  616.                NCVT = 0
  617.                IF( WNTVO .OR. WNTVAS ) THEN
  618. *
  619. *                 If right singular vectors desired, generate P'.
  620. *                 (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
  621. *                 (RWorkspace: 0)
  622. *
  623.                   CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
  624.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  625.                   NCVT = N
  626.                END IF
  627.                IRWORK = IE + N
  628. *
  629. *              Perform bidiagonal QR iteration, computing right
  630. *              singular vectors of A in A if desired
  631. *              (CWorkspace: 0)
  632. *              (RWorkspace: need BDSPAC)
  633. *
  634.                CALL ZBDSQR( 'U', N, NCVT, 0, 0, S, RWORK( IE ), A, LDA,
  635.      $                      CDUM, 1, CDUM, 1, RWORK( IRWORK ), INFO )
  636. *
  637. *              If right singular vectors desired in VT, copy them there
  638. *
  639.                IF( WNTVAS )
  640.      $            CALL ZLACPY( 'F', N, N, A, LDA, VT, LDVT )
  641. *
  642.             ELSE IF( WNTUO .AND. WNTVN ) THEN
  643. *
  644. *              Path 2 (M much larger than N, JOBU='O', JOBVT='N')
  645. *              N left singular vectors to be overwritten on A and
  646. *              no right singular vectors to be computed
  647. *
  648.                IF( LWORK.GE.N*N+3*N ) THEN
  649. *
  650. *                 Sufficient workspace for a fast algorithm
  651. *
  652.                   IR = 1
  653.                   IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*N ) THEN
  654. *
  655. *                    WORK(IU) is LDA by N, WORK(IR) is LDA by N
  656. *
  657.                      LDWRKU = LDA
  658.                      LDWRKR = LDA
  659.                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+N*N ) THEN
  660. *
  661. *                    WORK(IU) is LDA by N, WORK(IR) is N by N
  662. *
  663.                      LDWRKU = LDA
  664.                      LDWRKR = N
  665.                   ELSE
  666. *
  667. *                    WORK(IU) is LDWRKU by N, WORK(IR) is N by N
  668. *
  669.                      LDWRKU = ( LWORK-N*N ) / N
  670.                      LDWRKR = N
  671.                   END IF
  672.                   ITAU = IR + LDWRKR*N
  673.                   IWORK = ITAU + N
  674. *
  675. *                 Compute A=Q*R
  676. *                 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
  677. *                 (RWorkspace: 0)
  678. *
  679.                   CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
  680.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  681. *
  682. *                 Copy R to WORK(IR) and zero out below it
  683. *
  684.                   CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
  685.                   CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
  686.      $                         WORK( IR+1 ), LDWRKR )
  687. *
  688. *                 Generate Q in A
  689. *                 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
  690. *                 (RWorkspace: 0)
  691. *
  692.                   CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
  693.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  694.                   IE = 1
  695.                   ITAUQ = ITAU
  696.                   ITAUP = ITAUQ + N
  697.                   IWORK = ITAUP + N
  698. *
  699. *                 Bidiagonalize R in WORK(IR)
  700. *                 (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
  701. *                 (RWorkspace: need N)
  702. *
  703.                   CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
  704.      $                         WORK( ITAUQ ), WORK( ITAUP ),
  705.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  706. *
  707. *                 Generate left vectors bidiagonalizing R
  708. *                 (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
  709. *                 (RWorkspace: need 0)
  710. *
  711.                   CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
  712.      $                         WORK( ITAUQ ), WORK( IWORK ),
  713.      $                         LWORK-IWORK+1, IERR )
  714.                   IRWORK = IE + N
  715. *
  716. *                 Perform bidiagonal QR iteration, computing left
  717. *                 singular vectors of R in WORK(IR)
  718. *                 (CWorkspace: need N*N)
  719. *                 (RWorkspace: need BDSPAC)
  720. *
  721.                   CALL ZBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM, 1,
  722.      $                         WORK( IR ), LDWRKR, CDUM, 1,
  723.      $                         RWORK( IRWORK ), INFO )
  724.                   IU = ITAUQ
  725. *
  726. *                 Multiply Q in A by left singular vectors of R in
  727. *                 WORK(IR), storing result in WORK(IU) and copying to A
  728. *                 (CWorkspace: need N*N+N, prefer N*N+M*N)
  729. *                 (RWorkspace: 0)
  730. *
  731.                   DO 10 I = 1, M, LDWRKU
  732.                      CHUNK = MIN( M-I+1, LDWRKU )
  733.                      CALL ZGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ),
  734.      $                           LDA, WORK( IR ), LDWRKR, CZERO,
  735.      $                           WORK( IU ), LDWRKU )
  736.                      CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
  737.      $                            A( I, 1 ), LDA )
  738.    10             CONTINUE
  739. *
  740.                ELSE
  741. *
  742. *                 Insufficient workspace for a fast algorithm
  743. *
  744.                   IE = 1
  745.                   ITAUQ = 1
  746.                   ITAUP = ITAUQ + N
  747.                   IWORK = ITAUP + N
  748. *
  749. *                 Bidiagonalize A
  750. *                 (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
  751. *                 (RWorkspace: N)
  752. *
  753.                   CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ),
  754.      $                         WORK( ITAUQ ), WORK( ITAUP ),
  755.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  756. *
  757. *                 Generate left vectors bidiagonalizing A
  758. *                 (CWorkspace: need 3*N, prefer 2*N+N*NB)
  759. *                 (RWorkspace: 0)
  760. *
  761.                   CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
  762.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  763.                   IRWORK = IE + N
  764. *
  765. *                 Perform bidiagonal QR iteration, computing left
  766. *                 singular vectors of A in A
  767. *                 (CWorkspace: need 0)
  768. *                 (RWorkspace: need BDSPAC)
  769. *
  770.                   CALL ZBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM, 1,
  771.      $                         A, LDA, CDUM, 1, RWORK( IRWORK ), INFO )
  772. *
  773.                END IF
  774. *
  775.             ELSE IF( WNTUO .AND. WNTVAS ) THEN
  776. *
  777. *              Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A')
  778. *              N left singular vectors to be overwritten on A and
  779. *              N right singular vectors to be computed in VT
  780. *
  781.                IF( LWORK.GE.N*N+3*N ) THEN
  782. *
  783. *                 Sufficient workspace for a fast algorithm
  784. *
  785.                   IR = 1
  786.                   IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*N ) THEN
  787. *
  788. *                    WORK(IU) is LDA by N and WORK(IR) is LDA by N
  789. *
  790.                      LDWRKU = LDA
  791.                      LDWRKR = LDA
  792.                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+N*N ) THEN
  793. *
  794. *                    WORK(IU) is LDA by N and WORK(IR) is N by N
  795. *
  796.                      LDWRKU = LDA
  797.                      LDWRKR = N
  798.                   ELSE
  799. *
  800. *                    WORK(IU) is LDWRKU by N and WORK(IR) is N by N
  801. *
  802.                      LDWRKU = ( LWORK-N*N ) / N
  803.                      LDWRKR = N
  804.                   END IF
  805.                   ITAU = IR + LDWRKR*N
  806.                   IWORK = ITAU + N
  807. *
  808. *                 Compute A=Q*R
  809. *                 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
  810. *                 (RWorkspace: 0)
  811. *
  812.                   CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
  813.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  814. *
  815. *                 Copy R to VT, zeroing out below it
  816. *
  817.                   CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
  818.                   CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, VT( 2, 1 ),
  819.      $                         LDVT )
  820. *
  821. *                 Generate Q in A
  822. *                 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
  823. *                 (RWorkspace: 0)
  824. *
  825.                   CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
  826.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  827.                   IE = 1
  828.                   ITAUQ = ITAU
  829.                   ITAUP = ITAUQ + N
  830.                   IWORK = ITAUP + N
  831. *
  832. *                 Bidiagonalize R in VT, copying result to WORK(IR)
  833. *                 (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
  834. *                 (RWorkspace: need N)
  835. *
  836.                   CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ),
  837.      $                         WORK( ITAUQ ), WORK( ITAUP ),
  838.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  839.                   CALL ZLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR )
  840. *
  841. *                 Generate left vectors bidiagonalizing R in WORK(IR)
  842. *                 (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
  843. *                 (RWorkspace: 0)
  844. *
  845.                   CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
  846.      $                         WORK( ITAUQ ), WORK( IWORK ),
  847.      $                         LWORK-IWORK+1, IERR )
  848. *
  849. *                 Generate right vectors bidiagonalizing R in VT
  850. *                 (CWorkspace: need N*N+3*N-1, prefer N*N+2*N+(N-1)*NB)
  851. *                 (RWorkspace: 0)
  852. *
  853.                   CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
  854.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  855.                   IRWORK = IE + N
  856. *
  857. *                 Perform bidiagonal QR iteration, computing left
  858. *                 singular vectors of R in WORK(IR) and computing right
  859. *                 singular vectors of R in VT
  860. *                 (CWorkspace: need N*N)
  861. *                 (RWorkspace: need BDSPAC)
  862. *
  863.                   CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT,
  864.      $                         LDVT, WORK( IR ), LDWRKR, CDUM, 1,
  865.      $                         RWORK( IRWORK ), INFO )
  866.                   IU = ITAUQ
  867. *
  868. *                 Multiply Q in A by left singular vectors of R in
  869. *                 WORK(IR), storing result in WORK(IU) and copying to A
  870. *                 (CWorkspace: need N*N+N, prefer N*N+M*N)
  871. *                 (RWorkspace: 0)
  872. *
  873.                   DO 20 I = 1, M, LDWRKU
  874.                      CHUNK = MIN( M-I+1, LDWRKU )
  875.                      CALL ZGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ),
  876.      $                           LDA, WORK( IR ), LDWRKR, CZERO,
  877.      $                           WORK( IU ), LDWRKU )
  878.                      CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
  879.      $                            A( I, 1 ), LDA )
  880.    20             CONTINUE
  881. *
  882.                ELSE
  883. *
  884. *                 Insufficient workspace for a fast algorithm
  885. *
  886.                   ITAU = 1
  887.                   IWORK = ITAU + N
  888. *
  889. *                 Compute A=Q*R
  890. *                 (CWorkspace: need 2*N, prefer N+N*NB)
  891. *                 (RWorkspace: 0)
  892. *
  893.                   CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
  894.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  895. *
  896. *                 Copy R to VT, zeroing out below it
  897. *
  898.                   CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
  899.                   CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, VT( 2, 1 ),
  900.      $                         LDVT )
  901. *
  902. *                 Generate Q in A
  903. *                 (CWorkspace: need 2*N, prefer N+N*NB)
  904. *                 (RWorkspace: 0)
  905. *
  906.                   CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
  907.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  908.                   IE = 1
  909.                   ITAUQ = ITAU
  910.                   ITAUP = ITAUQ + N
  911.                   IWORK = ITAUP + N
  912. *
  913. *                 Bidiagonalize R in VT
  914. *                 (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
  915. *                 (RWorkspace: N)
  916. *
  917.                   CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ),
  918.      $                         WORK( ITAUQ ), WORK( ITAUP ),
  919.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  920. *
  921. *                 Multiply Q in A by left vectors bidiagonalizing R
  922. *                 (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
  923. *                 (RWorkspace: 0)
  924. *
  925.                   CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
  926.      $                         WORK( ITAUQ ), A, LDA, WORK( IWORK ),
  927.      $                         LWORK-IWORK+1, IERR )
  928. *
  929. *                 Generate right vectors bidiagonalizing R in VT
  930. *                 (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
  931. *                 (RWorkspace: 0)
  932. *
  933.                   CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
  934.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  935.                   IRWORK = IE + N
  936. *
  937. *                 Perform bidiagonal QR iteration, computing left
  938. *                 singular vectors of A in A and computing right
  939. *                 singular vectors of A in VT
  940. *                 (CWorkspace: 0)
  941. *                 (RWorkspace: need BDSPAC)
  942. *
  943.                   CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT,
  944.      $                         LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ),
  945.      $                         INFO )
  946. *
  947.                END IF
  948. *
  949.             ELSE IF( WNTUS ) THEN
  950. *
  951.                IF( WNTVN ) THEN
  952. *
  953. *                 Path 4 (M much larger than N, JOBU='S', JOBVT='N')
  954. *                 N left singular vectors to be computed in U and
  955. *                 no right singular vectors to be computed
  956. *
  957.                   IF( LWORK.GE.N*N+3*N ) THEN
  958. *
  959. *                    Sufficient workspace for a fast algorithm
  960. *
  961.                      IR = 1
  962.                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
  963. *
  964. *                       WORK(IR) is LDA by N
  965. *
  966.                         LDWRKR = LDA
  967.                      ELSE
  968. *
  969. *                       WORK(IR) is N by N
  970. *
  971.                         LDWRKR = N
  972.                      END IF
  973.                      ITAU = IR + LDWRKR*N
  974.                      IWORK = ITAU + N
  975. *
  976. *                    Compute A=Q*R
  977. *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
  978. *                    (RWorkspace: 0)
  979. *
  980.                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
  981.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  982. *
  983. *                    Copy R to WORK(IR), zeroing out below it
  984. *
  985.                      CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ),
  986.      $                            LDWRKR )
  987.                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
  988.      $                            WORK( IR+1 ), LDWRKR )
  989. *
  990. *                    Generate Q in A
  991. *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
  992. *                    (RWorkspace: 0)
  993. *
  994.                      CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
  995.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  996.                      IE = 1
  997.                      ITAUQ = ITAU
  998.                      ITAUP = ITAUQ + N
  999.                      IWORK = ITAUP + N
  1000. *
  1001. *                    Bidiagonalize R in WORK(IR)
  1002. *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
  1003. *                    (RWorkspace: need N)
  1004. *
  1005.                      CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S,
  1006.      $                            RWORK( IE ), WORK( ITAUQ ),
  1007.      $                            WORK( ITAUP ), WORK( IWORK ),
  1008.      $                            LWORK-IWORK+1, IERR )
  1009. *
  1010. *                    Generate left vectors bidiagonalizing R in WORK(IR)
  1011. *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
  1012. *                    (RWorkspace: 0)
  1013. *
  1014.                      CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
  1015.      $                            WORK( ITAUQ ), WORK( IWORK ),
  1016.      $                            LWORK-IWORK+1, IERR )
  1017.                      IRWORK = IE + N
  1018. *
  1019. *                    Perform bidiagonal QR iteration, computing left
  1020. *                    singular vectors of R in WORK(IR)
  1021. *                    (CWorkspace: need N*N)
  1022. *                    (RWorkspace: need BDSPAC)
  1023. *
  1024.                      CALL ZBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM,
  1025.      $                            1, WORK( IR ), LDWRKR, CDUM, 1,
  1026.      $                            RWORK( IRWORK ), INFO )
  1027. *
  1028. *                    Multiply Q in A by left singular vectors of R in
  1029. *                    WORK(IR), storing result in U
  1030. *                    (CWorkspace: need N*N)
  1031. *                    (RWorkspace: 0)
  1032. *
  1033.                      CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA,
  1034.      $                           WORK( IR ), LDWRKR, CZERO, U, LDU )
  1035. *
  1036.                   ELSE
  1037. *
  1038. *                    Insufficient workspace for a fast algorithm
  1039. *
  1040.                      ITAU = 1
  1041.                      IWORK = ITAU + N
  1042. *
  1043. *                    Compute A=Q*R, copying result to U
  1044. *                    (CWorkspace: need 2*N, prefer N+N*NB)
  1045. *                    (RWorkspace: 0)
  1046. *
  1047.                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
  1048.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1049.                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
  1050. *
  1051. *                    Generate Q in U
  1052. *                    (CWorkspace: need 2*N, prefer N+N*NB)
  1053. *                    (RWorkspace: 0)
  1054. *
  1055.                      CALL ZUNGQR( M, N, N, U, LDU, WORK( ITAU ),
  1056.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1057.                      IE = 1
  1058.                      ITAUQ = ITAU
  1059.                      ITAUP = ITAUQ + N
  1060.                      IWORK = ITAUP + N
  1061. *
  1062. *                    Zero out below R in A
  1063. *
  1064.                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
  1065.      $                            A( 2, 1 ), LDA )
  1066. *
  1067. *                    Bidiagonalize R in A
  1068. *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
  1069. *                    (RWorkspace: need N)
  1070. *
  1071.                      CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ),
  1072.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  1073.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1074. *
  1075. *                    Multiply Q in U by left vectors bidiagonalizing R
  1076. *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
  1077. *                    (RWorkspace: 0)
  1078. *
  1079.                      CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
  1080.      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
  1081.      $                            LWORK-IWORK+1, IERR )
  1082.                      IRWORK = IE + N
  1083. *
  1084. *                    Perform bidiagonal QR iteration, computing left
  1085. *                    singular vectors of A in U
  1086. *                    (CWorkspace: 0)
  1087. *                    (RWorkspace: need BDSPAC)
  1088. *
  1089.                      CALL ZBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM,
  1090.      $                            1, U, LDU, CDUM, 1, RWORK( IRWORK ),
  1091.      $                            INFO )
  1092. *
  1093.                   END IF
  1094. *
  1095.                ELSE IF( WNTVO ) THEN
  1096. *
  1097. *                 Path 5 (M much larger than N, JOBU='S', JOBVT='O')
  1098. *                 N left singular vectors to be computed in U and
  1099. *                 N right singular vectors to be overwritten on A
  1100. *
  1101.                   IF( LWORK.GE.2*N*N+3*N ) THEN
  1102. *
  1103. *                    Sufficient workspace for a fast algorithm
  1104. *
  1105.                      IU = 1
  1106.                      IF( LWORK.GE.WRKBL+2*LDA*N ) THEN
  1107. *
  1108. *                       WORK(IU) is LDA by N and WORK(IR) is LDA by N
  1109. *
  1110.                         LDWRKU = LDA
  1111.                         IR = IU + LDWRKU*N
  1112.                         LDWRKR = LDA
  1113.                      ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN
  1114. *
  1115. *                       WORK(IU) is LDA by N and WORK(IR) is N by N
  1116. *
  1117.                         LDWRKU = LDA
  1118.                         IR = IU + LDWRKU*N
  1119.                         LDWRKR = N
  1120.                      ELSE
  1121. *
  1122. *                       WORK(IU) is N by N and WORK(IR) is N by N
  1123. *
  1124.                         LDWRKU = N
  1125.                         IR = IU + LDWRKU*N
  1126.                         LDWRKR = N
  1127.                      END IF
  1128.                      ITAU = IR + LDWRKR*N
  1129.                      IWORK = ITAU + N
  1130. *
  1131. *                    Compute A=Q*R
  1132. *                    (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
  1133. *                    (RWorkspace: 0)
  1134. *
  1135.                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
  1136.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1137. *
  1138. *                    Copy R to WORK(IU), zeroing out below it
  1139. *
  1140.                      CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ),
  1141.      $                            LDWRKU )
  1142.                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
  1143.      $                            WORK( IU+1 ), LDWRKU )
  1144. *
  1145. *                    Generate Q in A
  1146. *                    (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
  1147. *                    (RWorkspace: 0)
  1148. *
  1149.                      CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
  1150.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1151.                      IE = 1
  1152.                      ITAUQ = ITAU
  1153.                      ITAUP = ITAUQ + N
  1154.                      IWORK = ITAUP + N
  1155. *
  1156. *                    Bidiagonalize R in WORK(IU), copying result to
  1157. *                    WORK(IR)
  1158. *                    (CWorkspace: need   2*N*N+3*N,
  1159. *                                 prefer 2*N*N+2*N+2*N*NB)
  1160. *                    (RWorkspace: need   N)
  1161. *
  1162.                      CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S,
  1163.      $                            RWORK( IE ), WORK( ITAUQ ),
  1164.      $                            WORK( ITAUP ), WORK( IWORK ),
  1165.      $                            LWORK-IWORK+1, IERR )
  1166.                      CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU,
  1167.      $                            WORK( IR ), LDWRKR )
  1168. *
  1169. *                    Generate left bidiagonalizing vectors in WORK(IU)
  1170. *                    (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
  1171. *                    (RWorkspace: 0)
  1172. *
  1173.                      CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
  1174.      $                            WORK( ITAUQ ), WORK( IWORK ),
  1175.      $                            LWORK-IWORK+1, IERR )
  1176. *
  1177. *                    Generate right bidiagonalizing vectors in WORK(IR)
  1178. *                    (CWorkspace: need   2*N*N+3*N-1,
  1179. *                                 prefer 2*N*N+2*N+(N-1)*NB)
  1180. *                    (RWorkspace: 0)
  1181. *
  1182.                      CALL ZUNGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
  1183.      $                            WORK( ITAUP ), WORK( IWORK ),
  1184.      $                            LWORK-IWORK+1, IERR )
  1185.                      IRWORK = IE + N
  1186. *
  1187. *                    Perform bidiagonal QR iteration, computing left
  1188. *                    singular vectors of R in WORK(IU) and computing
  1189. *                    right singular vectors of R in WORK(IR)
  1190. *                    (CWorkspace: need 2*N*N)
  1191. *                    (RWorkspace: need BDSPAC)
  1192. *
  1193.                      CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ),
  1194.      $                            WORK( IR ), LDWRKR, WORK( IU ),
  1195.      $                            LDWRKU, CDUM, 1, RWORK( IRWORK ),
  1196.      $                            INFO )
  1197. *
  1198. *                    Multiply Q in A by left singular vectors of R in
  1199. *                    WORK(IU), storing result in U
  1200. *                    (CWorkspace: need N*N)
  1201. *                    (RWorkspace: 0)
  1202. *
  1203.                      CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA,
  1204.      $                           WORK( IU ), LDWRKU, CZERO, U, LDU )
  1205. *
  1206. *                    Copy right singular vectors of R to A
  1207. *                    (CWorkspace: need N*N)
  1208. *                    (RWorkspace: 0)
  1209. *
  1210.                      CALL ZLACPY( 'F', N, N, WORK( IR ), LDWRKR, A,
  1211.      $                            LDA )
  1212. *
  1213.                   ELSE
  1214. *
  1215. *                    Insufficient workspace for a fast algorithm
  1216. *
  1217.                      ITAU = 1
  1218.                      IWORK = ITAU + N
  1219. *
  1220. *                    Compute A=Q*R, copying result to U
  1221. *                    (CWorkspace: need 2*N, prefer N+N*NB)
  1222. *                    (RWorkspace: 0)
  1223. *
  1224.                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
  1225.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1226.                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
  1227. *
  1228. *                    Generate Q in U
  1229. *                    (CWorkspace: need 2*N, prefer N+N*NB)
  1230. *                    (RWorkspace: 0)
  1231. *
  1232.                      CALL ZUNGQR( M, N, N, U, LDU, WORK( ITAU ),
  1233.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1234.                      IE = 1
  1235.                      ITAUQ = ITAU
  1236.                      ITAUP = ITAUQ + N
  1237.                      IWORK = ITAUP + N
  1238. *
  1239. *                    Zero out below R in A
  1240. *
  1241.                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
  1242.      $                            A( 2, 1 ), LDA )
  1243. *
  1244. *                    Bidiagonalize R in A
  1245. *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
  1246. *                    (RWorkspace: need N)
  1247. *
  1248.                      CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ),
  1249.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  1250.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1251. *
  1252. *                    Multiply Q in U by left vectors bidiagonalizing R
  1253. *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
  1254. *                    (RWorkspace: 0)
  1255. *
  1256.                      CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
  1257.      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
  1258.      $                            LWORK-IWORK+1, IERR )
  1259. *
  1260. *                    Generate right vectors bidiagonalizing R in A
  1261. *                    (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
  1262. *                    (RWorkspace: 0)
  1263. *
  1264.                      CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
  1265.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1266.                      IRWORK = IE + N
  1267. *
  1268. *                    Perform bidiagonal QR iteration, computing left
  1269. *                    singular vectors of A in U and computing right
  1270. *                    singular vectors of A in A
  1271. *                    (CWorkspace: 0)
  1272. *                    (RWorkspace: need BDSPAC)
  1273. *
  1274.                      CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), A,
  1275.      $                            LDA, U, LDU, CDUM, 1, RWORK( IRWORK ),
  1276.      $                            INFO )
  1277. *
  1278.                   END IF
  1279. *
  1280.                ELSE IF( WNTVAS ) THEN
  1281. *
  1282. *                 Path 6 (M much larger than N, JOBU='S', JOBVT='S'
  1283. *                         or 'A')
  1284. *                 N left singular vectors to be computed in U and
  1285. *                 N right singular vectors to be computed in VT
  1286. *
  1287.                   IF( LWORK.GE.N*N+3*N ) THEN
  1288. *
  1289. *                    Sufficient workspace for a fast algorithm
  1290. *
  1291.                      IU = 1
  1292.                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
  1293. *
  1294. *                       WORK(IU) is LDA by N
  1295. *
  1296.                         LDWRKU = LDA
  1297.                      ELSE
  1298. *
  1299. *                       WORK(IU) is N by N
  1300. *
  1301.                         LDWRKU = N
  1302.                      END IF
  1303.                      ITAU = IU + LDWRKU*N
  1304.                      IWORK = ITAU + N
  1305. *
  1306. *                    Compute A=Q*R
  1307. *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
  1308. *                    (RWorkspace: 0)
  1309. *
  1310.                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
  1311.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1312. *
  1313. *                    Copy R to WORK(IU), zeroing out below it
  1314. *
  1315.                      CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ),
  1316.      $                            LDWRKU )
  1317.                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
  1318.      $                            WORK( IU+1 ), LDWRKU )
  1319. *
  1320. *                    Generate Q in A
  1321. *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
  1322. *                    (RWorkspace: 0)
  1323. *
  1324.                      CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
  1325.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1326.                      IE = 1
  1327.                      ITAUQ = ITAU
  1328.                      ITAUP = ITAUQ + N
  1329.                      IWORK = ITAUP + N
  1330. *
  1331. *                    Bidiagonalize R in WORK(IU), copying result to VT
  1332. *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
  1333. *                    (RWorkspace: need N)
  1334. *
  1335.                      CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S,
  1336.      $                            RWORK( IE ), WORK( ITAUQ ),
  1337.      $                            WORK( ITAUP ), WORK( IWORK ),
  1338.      $                            LWORK-IWORK+1, IERR )
  1339.                      CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT,
  1340.      $                            LDVT )
  1341. *
  1342. *                    Generate left bidiagonalizing vectors in WORK(IU)
  1343. *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
  1344. *                    (RWorkspace: 0)
  1345. *
  1346.                      CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
  1347.      $                            WORK( ITAUQ ), WORK( IWORK ),
  1348.      $                            LWORK-IWORK+1, IERR )
  1349. *
  1350. *                    Generate right bidiagonalizing vectors in VT
  1351. *                    (CWorkspace: need   N*N+3*N-1,
  1352. *                                 prefer N*N+2*N+(N-1)*NB)
  1353. *                    (RWorkspace: 0)
  1354. *
  1355.                      CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
  1356.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1357.                      IRWORK = IE + N
  1358. *
  1359. *                    Perform bidiagonal QR iteration, computing left
  1360. *                    singular vectors of R in WORK(IU) and computing
  1361. *                    right singular vectors of R in VT
  1362. *                    (CWorkspace: need N*N)
  1363. *                    (RWorkspace: need BDSPAC)
  1364. *
  1365.                      CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT,
  1366.      $                            LDVT, WORK( IU ), LDWRKU, CDUM, 1,
  1367.      $                            RWORK( IRWORK ), INFO )
  1368. *
  1369. *                    Multiply Q in A by left singular vectors of R in
  1370. *                    WORK(IU), storing result in U
  1371. *                    (CWorkspace: need N*N)
  1372. *                    (RWorkspace: 0)
  1373. *
  1374.                      CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA,
  1375.      $                           WORK( IU ), LDWRKU, CZERO, U, LDU )
  1376. *
  1377.                   ELSE
  1378. *
  1379. *                    Insufficient workspace for a fast algorithm
  1380. *
  1381.                      ITAU = 1
  1382.                      IWORK = ITAU + N
  1383. *
  1384. *                    Compute A=Q*R, copying result to U
  1385. *                    (CWorkspace: need 2*N, prefer N+N*NB)
  1386. *                    (RWorkspace: 0)
  1387. *
  1388.                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
  1389.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1390.                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
  1391. *
  1392. *                    Generate Q in U
  1393. *                    (CWorkspace: need 2*N, prefer N+N*NB)
  1394. *                    (RWorkspace: 0)
  1395. *
  1396.                      CALL ZUNGQR( M, N, N, U, LDU, WORK( ITAU ),
  1397.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1398. *
  1399. *                    Copy R to VT, zeroing out below it
  1400. *
  1401.                      CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
  1402.                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
  1403.      $                            VT( 2, 1 ), LDVT )
  1404.                      IE = 1
  1405.                      ITAUQ = ITAU
  1406.                      ITAUP = ITAUQ + N
  1407.                      IWORK = ITAUP + N
  1408. *
  1409. *                    Bidiagonalize R in VT
  1410. *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
  1411. *                    (RWorkspace: need N)
  1412. *
  1413.                      CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ),
  1414.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  1415.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1416. *
  1417. *                    Multiply Q in U by left bidiagonalizing vectors
  1418. *                    in VT
  1419. *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
  1420. *                    (RWorkspace: 0)
  1421. *
  1422.                      CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
  1423.      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
  1424.      $                            LWORK-IWORK+1, IERR )
  1425. *
  1426. *                    Generate right bidiagonalizing vectors in VT
  1427. *                    (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
  1428. *                    (RWorkspace: 0)
  1429. *
  1430.                      CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
  1431.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1432.                      IRWORK = IE + N
  1433. *
  1434. *                    Perform bidiagonal QR iteration, computing left
  1435. *                    singular vectors of A in U and computing right
  1436. *                    singular vectors of A in VT
  1437. *                    (CWorkspace: 0)
  1438. *                    (RWorkspace: need BDSPAC)
  1439. *
  1440.                      CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT,
  1441.      $                            LDVT, U, LDU, CDUM, 1,
  1442.      $                            RWORK( IRWORK ), INFO )
  1443. *
  1444.                   END IF
  1445. *
  1446.                END IF
  1447. *
  1448.             ELSE IF( WNTUA ) THEN
  1449. *
  1450.                IF( WNTVN ) THEN
  1451. *
  1452. *                 Path 7 (M much larger than N, JOBU='A', JOBVT='N')
  1453. *                 M left singular vectors to be computed in U and
  1454. *                 no right singular vectors to be computed
  1455. *
  1456.                   IF( LWORK.GE.N*N+MAX( N+M, 3*N ) ) THEN
  1457. *
  1458. *                    Sufficient workspace for a fast algorithm
  1459. *
  1460.                      IR = 1
  1461.                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
  1462. *
  1463. *                       WORK(IR) is LDA by N
  1464. *
  1465.                         LDWRKR = LDA
  1466.                      ELSE
  1467. *
  1468. *                       WORK(IR) is N by N
  1469. *
  1470.                         LDWRKR = N
  1471.                      END IF
  1472.                      ITAU = IR + LDWRKR*N
  1473.                      IWORK = ITAU + N
  1474. *
  1475. *                    Compute A=Q*R, copying result to U
  1476. *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
  1477. *                    (RWorkspace: 0)
  1478. *
  1479.                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
  1480.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1481.                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
  1482. *
  1483. *                    Copy R to WORK(IR), zeroing out below it
  1484. *
  1485.                      CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ),
  1486.      $                            LDWRKR )
  1487.                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
  1488.      $                            WORK( IR+1 ), LDWRKR )
  1489. *
  1490. *                    Generate Q in U
  1491. *                    (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB)
  1492. *                    (RWorkspace: 0)
  1493. *
  1494.                      CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
  1495.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1496.                      IE = 1
  1497.                      ITAUQ = ITAU
  1498.                      ITAUP = ITAUQ + N
  1499.                      IWORK = ITAUP + N
  1500. *
  1501. *                    Bidiagonalize R in WORK(IR)
  1502. *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
  1503. *                    (RWorkspace: need N)
  1504. *
  1505.                      CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S,
  1506.      $                            RWORK( IE ), WORK( ITAUQ ),
  1507.      $                            WORK( ITAUP ), WORK( IWORK ),
  1508.      $                            LWORK-IWORK+1, IERR )
  1509. *
  1510. *                    Generate left bidiagonalizing vectors in WORK(IR)
  1511. *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
  1512. *                    (RWorkspace: 0)
  1513. *
  1514.                      CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
  1515.      $                            WORK( ITAUQ ), WORK( IWORK ),
  1516.      $                            LWORK-IWORK+1, IERR )
  1517.                      IRWORK = IE + N
  1518. *
  1519. *                    Perform bidiagonal QR iteration, computing left
  1520. *                    singular vectors of R in WORK(IR)
  1521. *                    (CWorkspace: need N*N)
  1522. *                    (RWorkspace: need BDSPAC)
  1523. *
  1524.                      CALL ZBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM,
  1525.      $                            1, WORK( IR ), LDWRKR, CDUM, 1,
  1526.      $                            RWORK( IRWORK ), INFO )
  1527. *
  1528. *                    Multiply Q in U by left singular vectors of R in
  1529. *                    WORK(IR), storing result in A
  1530. *                    (CWorkspace: need N*N)
  1531. *                    (RWorkspace: 0)
  1532. *
  1533.                      CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU,
  1534.      $                           WORK( IR ), LDWRKR, CZERO, A, LDA )
  1535. *
  1536. *                    Copy left singular vectors of A from A to U
  1537. *
  1538.                      CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
  1539. *
  1540.                   ELSE
  1541. *
  1542. *                    Insufficient workspace for a fast algorithm
  1543. *
  1544.                      ITAU = 1
  1545.                      IWORK = ITAU + N
  1546. *
  1547. *                    Compute A=Q*R, copying result to U
  1548. *                    (CWorkspace: need 2*N, prefer N+N*NB)
  1549. *                    (RWorkspace: 0)
  1550. *
  1551.                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
  1552.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1553.                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
  1554. *
  1555. *                    Generate Q in U
  1556. *                    (CWorkspace: need N+M, prefer N+M*NB)
  1557. *                    (RWorkspace: 0)
  1558. *
  1559.                      CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
  1560.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1561.                      IE = 1
  1562.                      ITAUQ = ITAU
  1563.                      ITAUP = ITAUQ + N
  1564.                      IWORK = ITAUP + N
  1565. *
  1566. *                    Zero out below R in A
  1567. *
  1568.                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
  1569.      $                            A( 2, 1 ), LDA )
  1570. *
  1571. *                    Bidiagonalize R in A
  1572. *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
  1573. *                    (RWorkspace: need N)
  1574. *
  1575.                      CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ),
  1576.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  1577.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1578. *
  1579. *                    Multiply Q in U by left bidiagonalizing vectors
  1580. *                    in A
  1581. *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
  1582. *                    (RWorkspace: 0)
  1583. *
  1584.                      CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
  1585.      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
  1586.      $                            LWORK-IWORK+1, IERR )
  1587.                      IRWORK = IE + N
  1588. *
  1589. *                    Perform bidiagonal QR iteration, computing left
  1590. *                    singular vectors of A in U
  1591. *                    (CWorkspace: 0)
  1592. *                    (RWorkspace: need BDSPAC)
  1593. *
  1594.                      CALL ZBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM,
  1595.      $                            1, U, LDU, CDUM, 1, RWORK( IRWORK ),
  1596.      $                            INFO )
  1597. *
  1598.                   END IF
  1599. *
  1600.                ELSE IF( WNTVO ) THEN
  1601. *
  1602. *                 Path 8 (M much larger than N, JOBU='A', JOBVT='O')
  1603. *                 M left singular vectors to be computed in U and
  1604. *                 N right singular vectors to be overwritten on A
  1605. *
  1606.                   IF( LWORK.GE.2*N*N+MAX( N+M, 3*N ) ) THEN
  1607. *
  1608. *                    Sufficient workspace for a fast algorithm
  1609. *
  1610.                      IU = 1
  1611.                      IF( LWORK.GE.WRKBL+2*LDA*N ) THEN
  1612. *
  1613. *                       WORK(IU) is LDA by N and WORK(IR) is LDA by N
  1614. *
  1615.                         LDWRKU = LDA
  1616.                         IR = IU + LDWRKU*N
  1617.                         LDWRKR = LDA
  1618.                      ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN
  1619. *
  1620. *                       WORK(IU) is LDA by N and WORK(IR) is N by N
  1621. *
  1622.                         LDWRKU = LDA
  1623.                         IR = IU + LDWRKU*N
  1624.                         LDWRKR = N
  1625.                      ELSE
  1626. *
  1627. *                       WORK(IU) is N by N and WORK(IR) is N by N
  1628. *
  1629.                         LDWRKU = N
  1630.                         IR = IU + LDWRKU*N
  1631.                         LDWRKR = N
  1632.                      END IF
  1633.                      ITAU = IR + LDWRKR*N
  1634.                      IWORK = ITAU + N
  1635. *
  1636. *                    Compute A=Q*R, copying result to U
  1637. *                    (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
  1638. *                    (RWorkspace: 0)
  1639. *
  1640.                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
  1641.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1642.                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
  1643. *
  1644. *                    Generate Q in U
  1645. *                    (CWorkspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB)
  1646. *                    (RWorkspace: 0)
  1647. *
  1648.                      CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
  1649.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1650. *
  1651. *                    Copy R to WORK(IU), zeroing out below it
  1652. *
  1653.                      CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ),
  1654.      $                            LDWRKU )
  1655.                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
  1656.      $                            WORK( IU+1 ), LDWRKU )
  1657.                      IE = 1
  1658.                      ITAUQ = ITAU
  1659.                      ITAUP = ITAUQ + N
  1660.                      IWORK = ITAUP + N
  1661. *
  1662. *                    Bidiagonalize R in WORK(IU), copying result to
  1663. *                    WORK(IR)
  1664. *                    (CWorkspace: need   2*N*N+3*N,
  1665. *                                 prefer 2*N*N+2*N+2*N*NB)
  1666. *                    (RWorkspace: need   N)
  1667. *
  1668.                      CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S,
  1669.      $                            RWORK( IE ), WORK( ITAUQ ),
  1670.      $                            WORK( ITAUP ), WORK( IWORK ),
  1671.      $                            LWORK-IWORK+1, IERR )
  1672.                      CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU,
  1673.      $                            WORK( IR ), LDWRKR )
  1674. *
  1675. *                    Generate left bidiagonalizing vectors in WORK(IU)
  1676. *                    (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
  1677. *                    (RWorkspace: 0)
  1678. *
  1679.                      CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
  1680.      $                            WORK( ITAUQ ), WORK( IWORK ),
  1681.      $                            LWORK-IWORK+1, IERR )
  1682. *
  1683. *                    Generate right bidiagonalizing vectors in WORK(IR)
  1684. *                    (CWorkspace: need   2*N*N+3*N-1,
  1685. *                                 prefer 2*N*N+2*N+(N-1)*NB)
  1686. *                    (RWorkspace: 0)
  1687. *
  1688.                      CALL ZUNGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
  1689.      $                            WORK( ITAUP ), WORK( IWORK ),
  1690.      $                            LWORK-IWORK+1, IERR )
  1691.                      IRWORK = IE + N
  1692. *
  1693. *                    Perform bidiagonal QR iteration, computing left
  1694. *                    singular vectors of R in WORK(IU) and computing
  1695. *                    right singular vectors of R in WORK(IR)
  1696. *                    (CWorkspace: need 2*N*N)
  1697. *                    (RWorkspace: need BDSPAC)
  1698. *
  1699.                      CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ),
  1700.      $                            WORK( IR ), LDWRKR, WORK( IU ),
  1701.      $                            LDWRKU, CDUM, 1, RWORK( IRWORK ),
  1702.      $                            INFO )
  1703. *
  1704. *                    Multiply Q in U by left singular vectors of R in
  1705. *                    WORK(IU), storing result in A
  1706. *                    (CWorkspace: need N*N)
  1707. *                    (RWorkspace: 0)
  1708. *
  1709.                      CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU,
  1710.      $                           WORK( IU ), LDWRKU, CZERO, A, LDA )
  1711. *
  1712. *                    Copy left singular vectors of A from A to U
  1713. *
  1714.                      CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
  1715. *
  1716. *                    Copy right singular vectors of R from WORK(IR) to A
  1717. *
  1718.                      CALL ZLACPY( 'F', N, N, WORK( IR ), LDWRKR, A,
  1719.      $                            LDA )
  1720. *
  1721.                   ELSE
  1722. *
  1723. *                    Insufficient workspace for a fast algorithm
  1724. *
  1725.                      ITAU = 1
  1726.                      IWORK = ITAU + N
  1727. *
  1728. *                    Compute A=Q*R, copying result to U
  1729. *                    (CWorkspace: need 2*N, prefer N+N*NB)
  1730. *                    (RWorkspace: 0)
  1731. *
  1732.                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
  1733.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1734.                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
  1735. *
  1736. *                    Generate Q in U
  1737. *                    (CWorkspace: need N+M, prefer N+M*NB)
  1738. *                    (RWorkspace: 0)
  1739. *
  1740.                      CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
  1741.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1742.                      IE = 1
  1743.                      ITAUQ = ITAU
  1744.                      ITAUP = ITAUQ + N
  1745.                      IWORK = ITAUP + N
  1746. *
  1747. *                    Zero out below R in A
  1748. *
  1749.                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
  1750.      $                            A( 2, 1 ), LDA )
  1751. *
  1752. *                    Bidiagonalize R in A
  1753. *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
  1754. *                    (RWorkspace: need N)
  1755. *
  1756.                      CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ),
  1757.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  1758.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1759. *
  1760. *                    Multiply Q in U by left bidiagonalizing vectors
  1761. *                    in A
  1762. *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
  1763. *                    (RWorkspace: 0)
  1764. *
  1765.                      CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
  1766.      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
  1767.      $                            LWORK-IWORK+1, IERR )
  1768. *
  1769. *                    Generate right bidiagonalizing vectors in A
  1770. *                    (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
  1771. *                    (RWorkspace: 0)
  1772. *
  1773.                      CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
  1774.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1775.                      IRWORK = IE + N
  1776. *
  1777. *                    Perform bidiagonal QR iteration, computing left
  1778. *                    singular vectors of A in U and computing right
  1779. *                    singular vectors of A in A
  1780. *                    (CWorkspace: 0)
  1781. *                    (RWorkspace: need BDSPAC)
  1782. *
  1783.                      CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), A,
  1784.      $                            LDA, U, LDU, CDUM, 1, RWORK( IRWORK ),
  1785.      $                            INFO )
  1786. *
  1787.                   END IF
  1788. *
  1789.                ELSE IF( WNTVAS ) THEN
  1790. *
  1791. *                 Path 9 (M much larger than N, JOBU='A', JOBVT='S'
  1792. *                         or 'A')
  1793. *                 M left singular vectors to be computed in U and
  1794. *                 N right singular vectors to be computed in VT
  1795. *
  1796.                   IF( LWORK.GE.N*N+MAX( N+M, 3*N ) ) THEN
  1797. *
  1798. *                    Sufficient workspace for a fast algorithm
  1799. *
  1800.                      IU = 1
  1801.                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
  1802. *
  1803. *                       WORK(IU) is LDA by N
  1804. *
  1805.                         LDWRKU = LDA
  1806.                      ELSE
  1807. *
  1808. *                       WORK(IU) is N by N
  1809. *
  1810.                         LDWRKU = N
  1811.                      END IF
  1812.                      ITAU = IU + LDWRKU*N
  1813.                      IWORK = ITAU + N
  1814. *
  1815. *                    Compute A=Q*R, copying result to U
  1816. *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
  1817. *                    (RWorkspace: 0)
  1818. *
  1819.                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
  1820.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1821.                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
  1822. *
  1823. *                    Generate Q in U
  1824. *                    (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB)
  1825. *                    (RWorkspace: 0)
  1826. *
  1827.                      CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
  1828.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1829. *
  1830. *                    Copy R to WORK(IU), zeroing out below it
  1831. *
  1832.                      CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ),
  1833.      $                            LDWRKU )
  1834.                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
  1835.      $                            WORK( IU+1 ), LDWRKU )
  1836.                      IE = 1
  1837.                      ITAUQ = ITAU
  1838.                      ITAUP = ITAUQ + N
  1839.                      IWORK = ITAUP + N
  1840. *
  1841. *                    Bidiagonalize R in WORK(IU), copying result to VT
  1842. *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
  1843. *                    (RWorkspace: need N)
  1844. *
  1845.                      CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S,
  1846.      $                            RWORK( IE ), WORK( ITAUQ ),
  1847.      $                            WORK( ITAUP ), WORK( IWORK ),
  1848.      $                            LWORK-IWORK+1, IERR )
  1849.                      CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT,
  1850.      $                            LDVT )
  1851. *
  1852. *                    Generate left bidiagonalizing vectors in WORK(IU)
  1853. *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
  1854. *                    (RWorkspace: 0)
  1855. *
  1856.                      CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
  1857.      $                            WORK( ITAUQ ), WORK( IWORK ),
  1858.      $                            LWORK-IWORK+1, IERR )
  1859. *
  1860. *                    Generate right bidiagonalizing vectors in VT
  1861. *                    (CWorkspace: need   N*N+3*N-1,
  1862. *                                 prefer N*N+2*N+(N-1)*NB)
  1863. *                    (RWorkspace: need   0)
  1864. *
  1865.                      CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
  1866.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1867.                      IRWORK = IE + N
  1868. *
  1869. *                    Perform bidiagonal QR iteration, computing left
  1870. *                    singular vectors of R in WORK(IU) and computing
  1871. *                    right singular vectors of R in VT
  1872. *                    (CWorkspace: need N*N)
  1873. *                    (RWorkspace: need BDSPAC)
  1874. *
  1875.                      CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT,
  1876.      $                            LDVT, WORK( IU ), LDWRKU, CDUM, 1,
  1877.      $                            RWORK( IRWORK ), INFO )
  1878. *
  1879. *                    Multiply Q in U by left singular vectors of R in
  1880. *                    WORK(IU), storing result in A
  1881. *                    (CWorkspace: need N*N)
  1882. *                    (RWorkspace: 0)
  1883. *
  1884.                      CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU,
  1885.      $                           WORK( IU ), LDWRKU, CZERO, A, LDA )
  1886. *
  1887. *                    Copy left singular vectors of A from A to U
  1888. *
  1889.                      CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
  1890. *
  1891.                   ELSE
  1892. *
  1893. *                    Insufficient workspace for a fast algorithm
  1894. *
  1895.                      ITAU = 1
  1896.                      IWORK = ITAU + N
  1897. *
  1898. *                    Compute A=Q*R, copying result to U
  1899. *                    (CWorkspace: need 2*N, prefer N+N*NB)
  1900. *                    (RWorkspace: 0)
  1901. *
  1902.                      CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
  1903.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1904.                      CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
  1905. *
  1906. *                    Generate Q in U
  1907. *                    (CWorkspace: need N+M, prefer N+M*NB)
  1908. *                    (RWorkspace: 0)
  1909. *
  1910.                      CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
  1911.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1912. *
  1913. *                    Copy R from A to VT, zeroing out below it
  1914. *
  1915.                      CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
  1916.                      CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
  1917.      $                            VT( 2, 1 ), LDVT )
  1918.                      IE = 1
  1919.                      ITAUQ = ITAU
  1920.                      ITAUP = ITAUQ + N
  1921.                      IWORK = ITAUP + N
  1922. *
  1923. *                    Bidiagonalize R in VT
  1924. *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
  1925. *                    (RWorkspace: need N)
  1926. *
  1927.                      CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ),
  1928.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  1929.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1930. *
  1931. *                    Multiply Q in U by left bidiagonalizing vectors
  1932. *                    in VT
  1933. *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
  1934. *                    (RWorkspace: 0)
  1935. *
  1936.                      CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
  1937.      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
  1938.      $                            LWORK-IWORK+1, IERR )
  1939. *
  1940. *                    Generate right bidiagonalizing vectors in VT
  1941. *                    (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
  1942. *                    (RWorkspace: 0)
  1943. *
  1944.                      CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
  1945.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1946.                      IRWORK = IE + N
  1947. *
  1948. *                    Perform bidiagonal QR iteration, computing left
  1949. *                    singular vectors of A in U and computing right
  1950. *                    singular vectors of A in VT
  1951. *                    (CWorkspace: 0)
  1952. *                    (RWorkspace: need BDSPAC)
  1953. *
  1954.                      CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT,
  1955.      $                            LDVT, U, LDU, CDUM, 1,
  1956.      $                            RWORK( IRWORK ), INFO )
  1957. *
  1958.                   END IF
  1959. *
  1960.                END IF
  1961. *
  1962.             END IF
  1963. *
  1964.          ELSE
  1965. *
  1966. *           M .LT. MNTHR
  1967. *
  1968. *           Path 10 (M at least N, but not much larger)
  1969. *           Reduce to bidiagonal form without QR decomposition
  1970. *
  1971.             IE = 1
  1972.             ITAUQ = 1
  1973.             ITAUP = ITAUQ + N
  1974.             IWORK = ITAUP + N
  1975. *
  1976. *           Bidiagonalize A
  1977. *           (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
  1978. *           (RWorkspace: need N)
  1979. *
  1980.             CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
  1981.      $                   WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
  1982.      $                   IERR )
  1983.             IF( WNTUAS ) THEN
  1984. *
  1985. *              If left singular vectors desired in U, copy result to U
  1986. *              and generate left bidiagonalizing vectors in U
  1987. *              (CWorkspace: need 2*N+NCU, prefer 2*N+NCU*NB)
  1988. *              (RWorkspace: 0)
  1989. *
  1990.                CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
  1991.                IF( WNTUS )
  1992.      $            NCU = N
  1993.                IF( WNTUA )
  1994.      $            NCU = M
  1995.                CALL ZUNGBR( 'Q', M, NCU, N, U, LDU, WORK( ITAUQ ),
  1996.      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
  1997.             END IF
  1998.             IF( WNTVAS ) THEN
  1999. *
  2000. *              If right singular vectors desired in VT, copy result to
  2001. *              VT and generate right bidiagonalizing vectors in VT
  2002. *              (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
  2003. *              (RWorkspace: 0)
  2004. *
  2005.                CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
  2006.                CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
  2007.      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
  2008.             END IF
  2009.             IF( WNTUO ) THEN
  2010. *
  2011. *              If left singular vectors desired in A, generate left
  2012. *              bidiagonalizing vectors in A
  2013. *              (CWorkspace: need 3*N, prefer 2*N+N*NB)
  2014. *              (RWorkspace: 0)
  2015. *
  2016.                CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
  2017.      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
  2018.             END IF
  2019.             IF( WNTVO ) THEN
  2020. *
  2021. *              If right singular vectors desired in A, generate right
  2022. *              bidiagonalizing vectors in A
  2023. *              (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
  2024. *              (RWorkspace: 0)
  2025. *
  2026.                CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
  2027.      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
  2028.             END IF
  2029.             IRWORK = IE + N
  2030.             IF( WNTUAS .OR. WNTUO )
  2031.      $         NRU = M
  2032.             IF( WNTUN )
  2033.      $         NRU = 0
  2034.             IF( WNTVAS .OR. WNTVO )
  2035.      $         NCVT = N
  2036.             IF( WNTVN )
  2037.      $         NCVT = 0
  2038.             IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN
  2039. *
  2040. *              Perform bidiagonal QR iteration, if desired, computing
  2041. *              left singular vectors in U and computing right singular
  2042. *              vectors in VT
  2043. *              (CWorkspace: 0)
  2044. *              (RWorkspace: need BDSPAC)
  2045. *
  2046.                CALL ZBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), VT,
  2047.      $                      LDVT, U, LDU, CDUM, 1, RWORK( IRWORK ),
  2048.      $                      INFO )
  2049.             ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN
  2050. *
  2051. *              Perform bidiagonal QR iteration, if desired, computing
  2052. *              left singular vectors in U and computing right singular
  2053. *              vectors in A
  2054. *              (CWorkspace: 0)
  2055. *              (RWorkspace: need BDSPAC)
  2056. *
  2057.                CALL ZBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), A,
  2058.      $                      LDA, U, LDU, CDUM, 1, RWORK( IRWORK ),
  2059.      $                      INFO )
  2060.             ELSE
  2061. *
  2062. *              Perform bidiagonal QR iteration, if desired, computing
  2063. *              left singular vectors in A and computing right singular
  2064. *              vectors in VT
  2065. *              (CWorkspace: 0)
  2066. *              (RWorkspace: need BDSPAC)
  2067. *
  2068.                CALL ZBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), VT,
  2069.      $                      LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ),
  2070.      $                      INFO )
  2071.             END IF
  2072. *
  2073.          END IF
  2074. *
  2075.       ELSE
  2076. *
  2077. *        A has more columns than rows. If A has sufficiently more
  2078. *        columns than rows, first reduce using the LQ decomposition (if
  2079. *        sufficient workspace available)
  2080. *
  2081.          IF( N.GE.MNTHR ) THEN
  2082. *
  2083.             IF( WNTVN ) THEN
  2084. *
  2085. *              Path 1t(N much larger than M, JOBVT='N')
  2086. *              No right singular vectors to be computed
  2087. *
  2088.                ITAU = 1
  2089.                IWORK = ITAU + M
  2090. *
  2091. *              Compute A=L*Q
  2092. *              (CWorkspace: need 2*M, prefer M+M*NB)
  2093. *              (RWorkspace: 0)
  2094. *
  2095.                CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
  2096.      $                      LWORK-IWORK+1, IERR )
  2097. *
  2098. *              Zero out above L
  2099. *
  2100.                CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
  2101.      $                      LDA )
  2102.                IE = 1
  2103.                ITAUQ = 1
  2104.                ITAUP = ITAUQ + M
  2105.                IWORK = ITAUP + M
  2106. *
  2107. *              Bidiagonalize L in A
  2108. *              (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
  2109. *              (RWorkspace: need M)
  2110. *
  2111.                CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
  2112.      $                      WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
  2113.      $                      IERR )
  2114.                IF( WNTUO .OR. WNTUAS ) THEN
  2115. *
  2116. *                 If left singular vectors desired, generate Q
  2117. *                 (CWorkspace: need 3*M, prefer 2*M+M*NB)
  2118. *                 (RWorkspace: 0)
  2119. *
  2120.                   CALL ZUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
  2121.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2122.                END IF
  2123.                IRWORK = IE + M
  2124.                NRU = 0
  2125.                IF( WNTUO .OR. WNTUAS )
  2126.      $            NRU = M
  2127. *
  2128. *              Perform bidiagonal QR iteration, computing left singular
  2129. *              vectors of A in A if desired
  2130. *              (CWorkspace: 0)
  2131. *              (RWorkspace: need BDSPAC)
  2132. *
  2133.                CALL ZBDSQR( 'U', M, 0, NRU, 0, S, RWORK( IE ), CDUM, 1,
  2134.      $                      A, LDA, CDUM, 1, RWORK( IRWORK ), INFO )
  2135. *
  2136. *              If left singular vectors desired in U, copy them there
  2137. *
  2138.                IF( WNTUAS )
  2139.      $            CALL ZLACPY( 'F', M, M, A, LDA, U, LDU )
  2140. *
  2141.             ELSE IF( WNTVO .AND. WNTUN ) THEN
  2142. *
  2143. *              Path 2t(N much larger than M, JOBU='N', JOBVT='O')
  2144. *              M right singular vectors to be overwritten on A and
  2145. *              no left singular vectors to be computed
  2146. *
  2147.                IF( LWORK.GE.M*M+3*M ) THEN
  2148. *
  2149. *                 Sufficient workspace for a fast algorithm
  2150. *
  2151.                   IR = 1
  2152.                   IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*M ) THEN
  2153. *
  2154. *                    WORK(IU) is LDA by N and WORK(IR) is LDA by M
  2155. *
  2156.                      LDWRKU = LDA
  2157.                      CHUNK = N
  2158.                      LDWRKR = LDA
  2159.                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+M*M ) THEN
  2160. *
  2161. *                    WORK(IU) is LDA by N and WORK(IR) is M by M
  2162. *
  2163.                      LDWRKU = LDA
  2164.                      CHUNK = N
  2165.                      LDWRKR = M
  2166.                   ELSE
  2167. *
  2168. *                    WORK(IU) is M by CHUNK and WORK(IR) is M by M
  2169. *
  2170.                      LDWRKU = M
  2171.                      CHUNK = ( LWORK-M*M ) / M
  2172.                      LDWRKR = M
  2173.                   END IF
  2174.                   ITAU = IR + LDWRKR*M
  2175.                   IWORK = ITAU + M
  2176. *
  2177. *                 Compute A=L*Q
  2178. *                 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
  2179. *                 (RWorkspace: 0)
  2180. *
  2181.                   CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
  2182.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2183. *
  2184. *                 Copy L to WORK(IR) and zero out above it
  2185. *
  2186.                   CALL ZLACPY( 'L', M, M, A, LDA, WORK( IR ), LDWRKR )
  2187.                   CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
  2188.      $                         WORK( IR+LDWRKR ), LDWRKR )
  2189. *
  2190. *                 Generate Q in A
  2191. *                 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
  2192. *                 (RWorkspace: 0)
  2193. *
  2194.                   CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
  2195.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2196.                   IE = 1
  2197.                   ITAUQ = ITAU
  2198.                   ITAUP = ITAUQ + M
  2199.                   IWORK = ITAUP + M
  2200. *
  2201. *                 Bidiagonalize L in WORK(IR)
  2202. *                 (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
  2203. *                 (RWorkspace: need M)
  2204. *
  2205.                   CALL ZGEBRD( M, M, WORK( IR ), LDWRKR, S, RWORK( IE ),
  2206.      $                         WORK( ITAUQ ), WORK( ITAUP ),
  2207.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2208. *
  2209. *                 Generate right vectors bidiagonalizing L
  2210. *                 (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB)
  2211. *                 (RWorkspace: 0)
  2212. *
  2213.                   CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
  2214.      $                         WORK( ITAUP ), WORK( IWORK ),
  2215.      $                         LWORK-IWORK+1, IERR )
  2216.                   IRWORK = IE + M
  2217. *
  2218. *                 Perform bidiagonal QR iteration, computing right
  2219. *                 singular vectors of L in WORK(IR)
  2220. *                 (CWorkspace: need M*M)
  2221. *                 (RWorkspace: need BDSPAC)
  2222. *
  2223.                   CALL ZBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ),
  2224.      $                         WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1,
  2225.      $                         RWORK( IRWORK ), INFO )
  2226.                   IU = ITAUQ
  2227. *
  2228. *                 Multiply right singular vectors of L in WORK(IR) by Q
  2229. *                 in A, storing result in WORK(IU) and copying to A
  2230. *                 (CWorkspace: need M*M+M, prefer M*M+M*N)
  2231. *                 (RWorkspace: 0)
  2232. *
  2233.                   DO 30 I = 1, N, CHUNK
  2234.                      BLK = MIN( N-I+1, CHUNK )
  2235.                      CALL ZGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IR ),
  2236.      $                           LDWRKR, A( 1, I ), LDA, CZERO,
  2237.      $                           WORK( IU ), LDWRKU )
  2238.                      CALL ZLACPY( 'F', M, BLK, WORK( IU ), LDWRKU,
  2239.      $                            A( 1, I ), LDA )
  2240.    30             CONTINUE
  2241. *
  2242.                ELSE
  2243. *
  2244. *                 Insufficient workspace for a fast algorithm
  2245. *
  2246.                   IE = 1
  2247.                   ITAUQ = 1
  2248.                   ITAUP = ITAUQ + M
  2249.                   IWORK = ITAUP + M
  2250. *
  2251. *                 Bidiagonalize A
  2252. *                 (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
  2253. *                 (RWorkspace: need M)
  2254. *
  2255.                   CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ),
  2256.      $                         WORK( ITAUQ ), WORK( ITAUP ),
  2257.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2258. *
  2259. *                 Generate right vectors bidiagonalizing A
  2260. *                 (CWorkspace: need 3*M, prefer 2*M+M*NB)
  2261. *                 (RWorkspace: 0)
  2262. *
  2263.                   CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
  2264.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2265.                   IRWORK = IE + M
  2266. *
  2267. *                 Perform bidiagonal QR iteration, computing right
  2268. *                 singular vectors of A in A
  2269. *                 (CWorkspace: 0)
  2270. *                 (RWorkspace: need BDSPAC)
  2271. *
  2272.                   CALL ZBDSQR( 'L', M, N, 0, 0, S, RWORK( IE ), A, LDA,
  2273.      $                         CDUM, 1, CDUM, 1, RWORK( IRWORK ), INFO )
  2274. *
  2275.                END IF
  2276. *
  2277.             ELSE IF( WNTVO .AND. WNTUAS ) THEN
  2278. *
  2279. *              Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O')
  2280. *              M right singular vectors to be overwritten on A and
  2281. *              M left singular vectors to be computed in U
  2282. *
  2283.                IF( LWORK.GE.M*M+3*M ) THEN
  2284. *
  2285. *                 Sufficient workspace for a fast algorithm
  2286. *
  2287.                   IR = 1
  2288.                   IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*M ) THEN
  2289. *
  2290. *                    WORK(IU) is LDA by N and WORK(IR) is LDA by M
  2291. *
  2292.                      LDWRKU = LDA
  2293.                      CHUNK = N
  2294.                      LDWRKR = LDA
  2295.                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+M*M ) THEN
  2296. *
  2297. *                    WORK(IU) is LDA by N and WORK(IR) is M by M
  2298. *
  2299.                      LDWRKU = LDA
  2300.                      CHUNK = N
  2301.                      LDWRKR = M
  2302.                   ELSE
  2303. *
  2304. *                    WORK(IU) is M by CHUNK and WORK(IR) is M by M
  2305. *
  2306.                      LDWRKU = M
  2307.                      CHUNK = ( LWORK-M*M ) / M
  2308.                      LDWRKR = M
  2309.                   END IF
  2310.                   ITAU = IR + LDWRKR*M
  2311.                   IWORK = ITAU + M
  2312. *
  2313. *                 Compute A=L*Q
  2314. *                 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
  2315. *                 (RWorkspace: 0)
  2316. *
  2317.                   CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
  2318.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2319. *
  2320. *                 Copy L to U, zeroing about above it
  2321. *
  2322.                   CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
  2323.                   CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, U( 1, 2 ),
  2324.      $                         LDU )
  2325. *
  2326. *                 Generate Q in A
  2327. *                 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
  2328. *                 (RWorkspace: 0)
  2329. *
  2330.                   CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
  2331.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2332.                   IE = 1
  2333.                   ITAUQ = ITAU
  2334.                   ITAUP = ITAUQ + M
  2335.                   IWORK = ITAUP + M
  2336. *
  2337. *                 Bidiagonalize L in U, copying result to WORK(IR)
  2338. *                 (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
  2339. *                 (RWorkspace: need M)
  2340. *
  2341.                   CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ),
  2342.      $                         WORK( ITAUQ ), WORK( ITAUP ),
  2343.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2344.                   CALL ZLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR )
  2345. *
  2346. *                 Generate right vectors bidiagonalizing L in WORK(IR)
  2347. *                 (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB)
  2348. *                 (RWorkspace: 0)
  2349. *
  2350.                   CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
  2351.      $                         WORK( ITAUP ), WORK( IWORK ),
  2352.      $                         LWORK-IWORK+1, IERR )
  2353. *
  2354. *                 Generate left vectors bidiagonalizing L in U
  2355. *                 (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
  2356. *                 (RWorkspace: 0)
  2357. *
  2358.                   CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
  2359.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2360.                   IRWORK = IE + M
  2361. *
  2362. *                 Perform bidiagonal QR iteration, computing left
  2363. *                 singular vectors of L in U, and computing right
  2364. *                 singular vectors of L in WORK(IR)
  2365. *                 (CWorkspace: need M*M)
  2366. *                 (RWorkspace: need BDSPAC)
  2367. *
  2368.                   CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
  2369.      $                         WORK( IR ), LDWRKR, U, LDU, CDUM, 1,
  2370.      $                         RWORK( IRWORK ), INFO )
  2371.                   IU = ITAUQ
  2372. *
  2373. *                 Multiply right singular vectors of L in WORK(IR) by Q
  2374. *                 in A, storing result in WORK(IU) and copying to A
  2375. *                 (CWorkspace: need M*M+M, prefer M*M+M*N))
  2376. *                 (RWorkspace: 0)
  2377. *
  2378.                   DO 40 I = 1, N, CHUNK
  2379.                      BLK = MIN( N-I+1, CHUNK )
  2380.                      CALL ZGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IR ),
  2381.      $                           LDWRKR, A( 1, I ), LDA, CZERO,
  2382.      $                           WORK( IU ), LDWRKU )
  2383.                      CALL ZLACPY( 'F', M, BLK, WORK( IU ), LDWRKU,
  2384.      $                            A( 1, I ), LDA )
  2385.    40             CONTINUE
  2386. *
  2387.                ELSE
  2388. *
  2389. *                 Insufficient workspace for a fast algorithm
  2390. *
  2391.                   ITAU = 1
  2392.                   IWORK = ITAU + M
  2393. *
  2394. *                 Compute A=L*Q
  2395. *                 (CWorkspace: need 2*M, prefer M+M*NB)
  2396. *                 (RWorkspace: 0)
  2397. *
  2398.                   CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
  2399.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2400. *
  2401. *                 Copy L to U, zeroing out above it
  2402. *
  2403.                   CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
  2404.                   CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, U( 1, 2 ),
  2405.      $                         LDU )
  2406. *
  2407. *                 Generate Q in A
  2408. *                 (CWorkspace: need 2*M, prefer M+M*NB)
  2409. *                 (RWorkspace: 0)
  2410. *
  2411.                   CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
  2412.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2413.                   IE = 1
  2414.                   ITAUQ = ITAU
  2415.                   ITAUP = ITAUQ + M
  2416.                   IWORK = ITAUP + M
  2417. *
  2418. *                 Bidiagonalize L in U
  2419. *                 (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
  2420. *                 (RWorkspace: need M)
  2421. *
  2422.                   CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ),
  2423.      $                         WORK( ITAUQ ), WORK( ITAUP ),
  2424.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2425. *
  2426. *                 Multiply right vectors bidiagonalizing L by Q in A
  2427. *                 (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
  2428. *                 (RWorkspace: 0)
  2429. *
  2430.                   CALL ZUNMBR( 'P', 'L', 'C', M, N, M, U, LDU,
  2431.      $                         WORK( ITAUP ), A, LDA, WORK( IWORK ),
  2432.      $                         LWORK-IWORK+1, IERR )
  2433. *
  2434. *                 Generate left vectors bidiagonalizing L in U
  2435. *                 (CWorkspace: need 3*M, prefer 2*M+M*NB)
  2436. *                 (RWorkspace: 0)
  2437. *
  2438.                   CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
  2439.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2440.                   IRWORK = IE + M
  2441. *
  2442. *                 Perform bidiagonal QR iteration, computing left
  2443. *                 singular vectors of A in U and computing right
  2444. *                 singular vectors of A in A
  2445. *                 (CWorkspace: 0)
  2446. *                 (RWorkspace: need BDSPAC)
  2447. *
  2448.                   CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), A, LDA,
  2449.      $                         U, LDU, CDUM, 1, RWORK( IRWORK ), INFO )
  2450. *
  2451.                END IF
  2452. *
  2453.             ELSE IF( WNTVS ) THEN
  2454. *
  2455.                IF( WNTUN ) THEN
  2456. *
  2457. *                 Path 4t(N much larger than M, JOBU='N', JOBVT='S')
  2458. *                 M right singular vectors to be computed in VT and
  2459. *                 no left singular vectors to be computed
  2460. *
  2461.                   IF( LWORK.GE.M*M+3*M ) THEN
  2462. *
  2463. *                    Sufficient workspace for a fast algorithm
  2464. *
  2465.                      IR = 1
  2466.                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
  2467. *
  2468. *                       WORK(IR) is LDA by M
  2469. *
  2470.                         LDWRKR = LDA
  2471.                      ELSE
  2472. *
  2473. *                       WORK(IR) is M by M
  2474. *
  2475.                         LDWRKR = M
  2476.                      END IF
  2477.                      ITAU = IR + LDWRKR*M
  2478.                      IWORK = ITAU + M
  2479. *
  2480. *                    Compute A=L*Q
  2481. *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
  2482. *                    (RWorkspace: 0)
  2483. *
  2484.                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
  2485.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2486. *
  2487. *                    Copy L to WORK(IR), zeroing out above it
  2488. *
  2489.                      CALL ZLACPY( 'L', M, M, A, LDA, WORK( IR ),
  2490.      $                            LDWRKR )
  2491.                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
  2492.      $                            WORK( IR+LDWRKR ), LDWRKR )
  2493. *
  2494. *                    Generate Q in A
  2495. *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
  2496. *                    (RWorkspace: 0)
  2497. *
  2498.                      CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
  2499.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2500.                      IE = 1
  2501.                      ITAUQ = ITAU
  2502.                      ITAUP = ITAUQ + M
  2503.                      IWORK = ITAUP + M
  2504. *
  2505. *                    Bidiagonalize L in WORK(IR)
  2506. *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
  2507. *                    (RWorkspace: need M)
  2508. *
  2509.                      CALL ZGEBRD( M, M, WORK( IR ), LDWRKR, S,
  2510.      $                            RWORK( IE ), WORK( ITAUQ ),
  2511.      $                            WORK( ITAUP ), WORK( IWORK ),
  2512.      $                            LWORK-IWORK+1, IERR )
  2513. *
  2514. *                    Generate right vectors bidiagonalizing L in
  2515. *                    WORK(IR)
  2516. *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB)
  2517. *                    (RWorkspace: 0)
  2518. *
  2519.                      CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
  2520.      $                            WORK( ITAUP ), WORK( IWORK ),
  2521.      $                            LWORK-IWORK+1, IERR )
  2522.                      IRWORK = IE + M
  2523. *
  2524. *                    Perform bidiagonal QR iteration, computing right
  2525. *                    singular vectors of L in WORK(IR)
  2526. *                    (CWorkspace: need M*M)
  2527. *                    (RWorkspace: need BDSPAC)
  2528. *
  2529.                      CALL ZBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ),
  2530.      $                            WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1,
  2531.      $                            RWORK( IRWORK ), INFO )
  2532. *
  2533. *                    Multiply right singular vectors of L in WORK(IR) by
  2534. *                    Q in A, storing result in VT
  2535. *                    (CWorkspace: need M*M)
  2536. *                    (RWorkspace: 0)
  2537. *
  2538.                      CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IR ),
  2539.      $                           LDWRKR, A, LDA, CZERO, VT, LDVT )
  2540. *
  2541.                   ELSE
  2542. *
  2543. *                    Insufficient workspace for a fast algorithm
  2544. *
  2545.                      ITAU = 1
  2546.                      IWORK = ITAU + M
  2547. *
  2548. *                    Compute A=L*Q
  2549. *                    (CWorkspace: need 2*M, prefer M+M*NB)
  2550. *                    (RWorkspace: 0)
  2551. *
  2552.                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
  2553.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2554. *
  2555. *                    Copy result to VT
  2556. *
  2557.                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
  2558. *
  2559. *                    Generate Q in VT
  2560. *                    (CWorkspace: need 2*M, prefer M+M*NB)
  2561. *                    (RWorkspace: 0)
  2562. *
  2563.                      CALL ZUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
  2564.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2565.                      IE = 1
  2566.                      ITAUQ = ITAU
  2567.                      ITAUP = ITAUQ + M
  2568.                      IWORK = ITAUP + M
  2569. *
  2570. *                    Zero out above L in A
  2571. *
  2572.                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
  2573.      $                            A( 1, 2 ), LDA )
  2574. *
  2575. *                    Bidiagonalize L in A
  2576. *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
  2577. *                    (RWorkspace: need M)
  2578. *
  2579.                      CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ),
  2580.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  2581.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2582. *
  2583. *                    Multiply right vectors bidiagonalizing L by Q in VT
  2584. *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
  2585. *                    (RWorkspace: 0)
  2586. *
  2587.                      CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA,
  2588.      $                            WORK( ITAUP ), VT, LDVT,
  2589.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2590.                      IRWORK = IE + M
  2591. *
  2592. *                    Perform bidiagonal QR iteration, computing right
  2593. *                    singular vectors of A in VT
  2594. *                    (CWorkspace: 0)
  2595. *                    (RWorkspace: need BDSPAC)
  2596. *
  2597.                      CALL ZBDSQR( 'U', M, N, 0, 0, S, RWORK( IE ), VT,
  2598.      $                            LDVT, CDUM, 1, CDUM, 1,
  2599.      $                            RWORK( IRWORK ), INFO )
  2600. *
  2601.                   END IF
  2602. *
  2603.                ELSE IF( WNTUO ) THEN
  2604. *
  2605. *                 Path 5t(N much larger than M, JOBU='O', JOBVT='S')
  2606. *                 M right singular vectors to be computed in VT and
  2607. *                 M left singular vectors to be overwritten on A
  2608. *
  2609.                   IF( LWORK.GE.2*M*M+3*M ) THEN
  2610. *
  2611. *                    Sufficient workspace for a fast algorithm
  2612. *
  2613.                      IU = 1
  2614.                      IF( LWORK.GE.WRKBL+2*LDA*M ) THEN
  2615. *
  2616. *                       WORK(IU) is LDA by M and WORK(IR) is LDA by M
  2617. *
  2618.                         LDWRKU = LDA
  2619.                         IR = IU + LDWRKU*M
  2620.                         LDWRKR = LDA
  2621.                      ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN
  2622. *
  2623. *                       WORK(IU) is LDA by M and WORK(IR) is M by M
  2624. *
  2625.                         LDWRKU = LDA
  2626.                         IR = IU + LDWRKU*M
  2627.                         LDWRKR = M
  2628.                      ELSE
  2629. *
  2630. *                       WORK(IU) is M by M and WORK(IR) is M by M
  2631. *
  2632.                         LDWRKU = M
  2633.                         IR = IU + LDWRKU*M
  2634.                         LDWRKR = M
  2635.                      END IF
  2636.                      ITAU = IR + LDWRKR*M
  2637.                      IWORK = ITAU + M
  2638. *
  2639. *                    Compute A=L*Q
  2640. *                    (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
  2641. *                    (RWorkspace: 0)
  2642. *
  2643.                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
  2644.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2645. *
  2646. *                    Copy L to WORK(IU), zeroing out below it
  2647. *
  2648.                      CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ),
  2649.      $                            LDWRKU )
  2650.                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
  2651.      $                            WORK( IU+LDWRKU ), LDWRKU )
  2652. *
  2653. *                    Generate Q in A
  2654. *                    (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
  2655. *                    (RWorkspace: 0)
  2656. *
  2657.                      CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
  2658.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2659.                      IE = 1
  2660.                      ITAUQ = ITAU
  2661.                      ITAUP = ITAUQ + M
  2662.                      IWORK = ITAUP + M
  2663. *
  2664. *                    Bidiagonalize L in WORK(IU), copying result to
  2665. *                    WORK(IR)
  2666. *                    (CWorkspace: need   2*M*M+3*M,
  2667. *                                 prefer 2*M*M+2*M+2*M*NB)
  2668. *                    (RWorkspace: need   M)
  2669. *
  2670.                      CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S,
  2671.      $                            RWORK( IE ), WORK( ITAUQ ),
  2672.      $                            WORK( ITAUP ), WORK( IWORK ),
  2673.      $                            LWORK-IWORK+1, IERR )
  2674.                      CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU,
  2675.      $                            WORK( IR ), LDWRKR )
  2676. *
  2677. *                    Generate right bidiagonalizing vectors in WORK(IU)
  2678. *                    (CWorkspace: need   2*M*M+3*M-1,
  2679. *                                 prefer 2*M*M+2*M+(M-1)*NB)
  2680. *                    (RWorkspace: 0)
  2681. *
  2682.                      CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
  2683.      $                            WORK( ITAUP ), WORK( IWORK ),
  2684.      $                            LWORK-IWORK+1, IERR )
  2685. *
  2686. *                    Generate left bidiagonalizing vectors in WORK(IR)
  2687. *                    (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
  2688. *                    (RWorkspace: 0)
  2689. *
  2690.                      CALL ZUNGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
  2691.      $                            WORK( ITAUQ ), WORK( IWORK ),
  2692.      $                            LWORK-IWORK+1, IERR )
  2693.                      IRWORK = IE + M
  2694. *
  2695. *                    Perform bidiagonal QR iteration, computing left
  2696. *                    singular vectors of L in WORK(IR) and computing
  2697. *                    right singular vectors of L in WORK(IU)
  2698. *                    (CWorkspace: need 2*M*M)
  2699. *                    (RWorkspace: need BDSPAC)
  2700. *
  2701.                      CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
  2702.      $                            WORK( IU ), LDWRKU, WORK( IR ),
  2703.      $                            LDWRKR, CDUM, 1, RWORK( IRWORK ),
  2704.      $                            INFO )
  2705. *
  2706. *                    Multiply right singular vectors of L in WORK(IU) by
  2707. *                    Q in A, storing result in VT
  2708. *                    (CWorkspace: need M*M)
  2709. *                    (RWorkspace: 0)
  2710. *
  2711.                      CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ),
  2712.      $                           LDWRKU, A, LDA, CZERO, VT, LDVT )
  2713. *
  2714. *                    Copy left singular vectors of L to A
  2715. *                    (CWorkspace: need M*M)
  2716. *                    (RWorkspace: 0)
  2717. *
  2718.                      CALL ZLACPY( 'F', M, M, WORK( IR ), LDWRKR, A,
  2719.      $                            LDA )
  2720. *
  2721.                   ELSE
  2722. *
  2723. *                    Insufficient workspace for a fast algorithm
  2724. *
  2725.                      ITAU = 1
  2726.                      IWORK = ITAU + M
  2727. *
  2728. *                    Compute A=L*Q, copying result to VT
  2729. *                    (CWorkspace: need 2*M, prefer M+M*NB)
  2730. *                    (RWorkspace: 0)
  2731. *
  2732.                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
  2733.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2734.                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
  2735. *
  2736. *                    Generate Q in VT
  2737. *                    (CWorkspace: need 2*M, prefer M+M*NB)
  2738. *                    (RWorkspace: 0)
  2739. *
  2740.                      CALL ZUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
  2741.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2742.                      IE = 1
  2743.                      ITAUQ = ITAU
  2744.                      ITAUP = ITAUQ + M
  2745.                      IWORK = ITAUP + M
  2746. *
  2747. *                    Zero out above L in A
  2748. *
  2749.                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
  2750.      $                            A( 1, 2 ), LDA )
  2751. *
  2752. *                    Bidiagonalize L in A
  2753. *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
  2754. *                    (RWorkspace: need M)
  2755. *
  2756.                      CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ),
  2757.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  2758.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2759. *
  2760. *                    Multiply right vectors bidiagonalizing L by Q in VT
  2761. *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
  2762. *                    (RWorkspace: 0)
  2763. *
  2764.                      CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA,
  2765.      $                            WORK( ITAUP ), VT, LDVT,
  2766.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2767. *
  2768. *                    Generate left bidiagonalizing vectors of L in A
  2769. *                    (CWorkspace: need 3*M, prefer 2*M+M*NB)
  2770. *                    (RWorkspace: 0)
  2771. *
  2772.                      CALL ZUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
  2773.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2774.                      IRWORK = IE + M
  2775. *
  2776. *                    Perform bidiagonal QR iteration, computing left
  2777. *                    singular vectors of A in A and computing right
  2778. *                    singular vectors of A in VT
  2779. *                    (CWorkspace: 0)
  2780. *                    (RWorkspace: need BDSPAC)
  2781. *
  2782.                      CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT,
  2783.      $                            LDVT, A, LDA, CDUM, 1,
  2784.      $                            RWORK( IRWORK ), INFO )
  2785. *
  2786.                   END IF
  2787. *
  2788.                ELSE IF( WNTUAS ) THEN
  2789. *
  2790. *                 Path 6t(N much larger than M, JOBU='S' or 'A',
  2791. *                         JOBVT='S')
  2792. *                 M right singular vectors to be computed in VT and
  2793. *                 M left singular vectors to be computed in U
  2794. *
  2795.                   IF( LWORK.GE.M*M+3*M ) THEN
  2796. *
  2797. *                    Sufficient workspace for a fast algorithm
  2798. *
  2799.                      IU = 1
  2800.                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
  2801. *
  2802. *                       WORK(IU) is LDA by N
  2803. *
  2804.                         LDWRKU = LDA
  2805.                      ELSE
  2806. *
  2807. *                       WORK(IU) is LDA by M
  2808. *
  2809.                         LDWRKU = M
  2810.                      END IF
  2811.                      ITAU = IU + LDWRKU*M
  2812.                      IWORK = ITAU + M
  2813. *
  2814. *                    Compute A=L*Q
  2815. *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
  2816. *                    (RWorkspace: 0)
  2817. *
  2818.                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
  2819.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2820. *
  2821. *                    Copy L to WORK(IU), zeroing out above it
  2822. *
  2823.                      CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ),
  2824.      $                            LDWRKU )
  2825.                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
  2826.      $                            WORK( IU+LDWRKU ), LDWRKU )
  2827. *
  2828. *                    Generate Q in A
  2829. *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
  2830. *                    (RWorkspace: 0)
  2831. *
  2832.                      CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
  2833.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2834.                      IE = 1
  2835.                      ITAUQ = ITAU
  2836.                      ITAUP = ITAUQ + M
  2837.                      IWORK = ITAUP + M
  2838. *
  2839. *                    Bidiagonalize L in WORK(IU), copying result to U
  2840. *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
  2841. *                    (RWorkspace: need M)
  2842. *
  2843.                      CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S,
  2844.      $                            RWORK( IE ), WORK( ITAUQ ),
  2845.      $                            WORK( ITAUP ), WORK( IWORK ),
  2846.      $                            LWORK-IWORK+1, IERR )
  2847.                      CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU, U,
  2848.      $                            LDU )
  2849. *
  2850. *                    Generate right bidiagonalizing vectors in WORK(IU)
  2851. *                    (CWorkspace: need   M*M+3*M-1,
  2852. *                                 prefer M*M+2*M+(M-1)*NB)
  2853. *                    (RWorkspace: 0)
  2854. *
  2855.                      CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
  2856.      $                            WORK( ITAUP ), WORK( IWORK ),
  2857.      $                            LWORK-IWORK+1, IERR )
  2858. *
  2859. *                    Generate left bidiagonalizing vectors in U
  2860. *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
  2861. *                    (RWorkspace: 0)
  2862. *
  2863.                      CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
  2864.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2865.                      IRWORK = IE + M
  2866. *
  2867. *                    Perform bidiagonal QR iteration, computing left
  2868. *                    singular vectors of L in U and computing right
  2869. *                    singular vectors of L in WORK(IU)
  2870. *                    (CWorkspace: need M*M)
  2871. *                    (RWorkspace: need BDSPAC)
  2872. *
  2873.                      CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
  2874.      $                            WORK( IU ), LDWRKU, U, LDU, CDUM, 1,
  2875.      $                            RWORK( IRWORK ), INFO )
  2876. *
  2877. *                    Multiply right singular vectors of L in WORK(IU) by
  2878. *                    Q in A, storing result in VT
  2879. *                    (CWorkspace: need M*M)
  2880. *                    (RWorkspace: 0)
  2881. *
  2882.                      CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ),
  2883.      $                           LDWRKU, A, LDA, CZERO, VT, LDVT )
  2884. *
  2885.                   ELSE
  2886. *
  2887. *                    Insufficient workspace for a fast algorithm
  2888. *
  2889.                      ITAU = 1
  2890.                      IWORK = ITAU + M
  2891. *
  2892. *                    Compute A=L*Q, copying result to VT
  2893. *                    (CWorkspace: need 2*M, prefer M+M*NB)
  2894. *                    (RWorkspace: 0)
  2895. *
  2896.                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
  2897.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2898.                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
  2899. *
  2900. *                    Generate Q in VT
  2901. *                    (CWorkspace: need 2*M, prefer M+M*NB)
  2902. *                    (RWorkspace: 0)
  2903. *
  2904.                      CALL ZUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
  2905.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2906. *
  2907. *                    Copy L to U, zeroing out above it
  2908. *
  2909.                      CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
  2910.                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
  2911.      $                            U( 1, 2 ), LDU )
  2912.                      IE = 1
  2913.                      ITAUQ = ITAU
  2914.                      ITAUP = ITAUQ + M
  2915.                      IWORK = ITAUP + M
  2916. *
  2917. *                    Bidiagonalize L in U
  2918. *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
  2919. *                    (RWorkspace: need M)
  2920. *
  2921.                      CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ),
  2922.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  2923.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2924. *
  2925. *                    Multiply right bidiagonalizing vectors in U by Q
  2926. *                    in VT
  2927. *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
  2928. *                    (RWorkspace: 0)
  2929. *
  2930.                      CALL ZUNMBR( 'P', 'L', 'C', M, N, M, U, LDU,
  2931.      $                            WORK( ITAUP ), VT, LDVT,
  2932.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2933. *
  2934. *                    Generate left bidiagonalizing vectors in U
  2935. *                    (CWorkspace: need 3*M, prefer 2*M+M*NB)
  2936. *                    (RWorkspace: 0)
  2937. *
  2938.                      CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
  2939.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2940.                      IRWORK = IE + M
  2941. *
  2942. *                    Perform bidiagonal QR iteration, computing left
  2943. *                    singular vectors of A in U and computing right
  2944. *                    singular vectors of A in VT
  2945. *                    (CWorkspace: 0)
  2946. *                    (RWorkspace: need BDSPAC)
  2947. *
  2948.                      CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT,
  2949.      $                            LDVT, U, LDU, CDUM, 1,
  2950.      $                            RWORK( IRWORK ), INFO )
  2951. *
  2952.                   END IF
  2953. *
  2954.                END IF
  2955. *
  2956.             ELSE IF( WNTVA ) THEN
  2957. *
  2958.                IF( WNTUN ) THEN
  2959. *
  2960. *                 Path 7t(N much larger than M, JOBU='N', JOBVT='A')
  2961. *                 N right singular vectors to be computed in VT and
  2962. *                 no left singular vectors to be computed
  2963. *
  2964.                   IF( LWORK.GE.M*M+MAX( N+M, 3*M ) ) THEN
  2965. *
  2966. *                    Sufficient workspace for a fast algorithm
  2967. *
  2968.                      IR = 1
  2969.                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
  2970. *
  2971. *                       WORK(IR) is LDA by M
  2972. *
  2973.                         LDWRKR = LDA
  2974.                      ELSE
  2975. *
  2976. *                       WORK(IR) is M by M
  2977. *
  2978.                         LDWRKR = M
  2979.                      END IF
  2980.                      ITAU = IR + LDWRKR*M
  2981.                      IWORK = ITAU + M
  2982. *
  2983. *                    Compute A=L*Q, copying result to VT
  2984. *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
  2985. *                    (RWorkspace: 0)
  2986. *
  2987.                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
  2988.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2989.                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
  2990. *
  2991. *                    Copy L to WORK(IR), zeroing out above it
  2992. *
  2993.                      CALL ZLACPY( 'L', M, M, A, LDA, WORK( IR ),
  2994.      $                            LDWRKR )
  2995.                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
  2996.      $                            WORK( IR+LDWRKR ), LDWRKR )
  2997. *
  2998. *                    Generate Q in VT
  2999. *                    (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB)
  3000. *                    (RWorkspace: 0)
  3001. *
  3002.                      CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
  3003.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3004.                      IE = 1
  3005.                      ITAUQ = ITAU
  3006.                      ITAUP = ITAUQ + M
  3007.                      IWORK = ITAUP + M
  3008. *
  3009. *                    Bidiagonalize L in WORK(IR)
  3010. *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
  3011. *                    (RWorkspace: need M)
  3012. *
  3013.                      CALL ZGEBRD( M, M, WORK( IR ), LDWRKR, S,
  3014.      $                            RWORK( IE ), WORK( ITAUQ ),
  3015.      $                            WORK( ITAUP ), WORK( IWORK ),
  3016.      $                            LWORK-IWORK+1, IERR )
  3017. *
  3018. *                    Generate right bidiagonalizing vectors in WORK(IR)
  3019. *                    (CWorkspace: need   M*M+3*M-1,
  3020. *                                 prefer M*M+2*M+(M-1)*NB)
  3021. *                    (RWorkspace: 0)
  3022. *
  3023.                      CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
  3024.      $                            WORK( ITAUP ), WORK( IWORK ),
  3025.      $                            LWORK-IWORK+1, IERR )
  3026.                      IRWORK = IE + M
  3027. *
  3028. *                    Perform bidiagonal QR iteration, computing right
  3029. *                    singular vectors of L in WORK(IR)
  3030. *                    (CWorkspace: need M*M)
  3031. *                    (RWorkspace: need BDSPAC)
  3032. *
  3033.                      CALL ZBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ),
  3034.      $                            WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1,
  3035.      $                            RWORK( IRWORK ), INFO )
  3036. *
  3037. *                    Multiply right singular vectors of L in WORK(IR) by
  3038. *                    Q in VT, storing result in A
  3039. *                    (CWorkspace: need M*M)
  3040. *                    (RWorkspace: 0)
  3041. *
  3042.                      CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IR ),
  3043.      $                           LDWRKR, VT, LDVT, CZERO, A, LDA )
  3044. *
  3045. *                    Copy right singular vectors of A from A to VT
  3046. *
  3047.                      CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
  3048. *
  3049.                   ELSE
  3050. *
  3051. *                    Insufficient workspace for a fast algorithm
  3052. *
  3053.                      ITAU = 1
  3054.                      IWORK = ITAU + M
  3055. *
  3056. *                    Compute A=L*Q, copying result to VT
  3057. *                    (CWorkspace: need 2*M, prefer M+M*NB)
  3058. *                    (RWorkspace: 0)
  3059. *
  3060.                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
  3061.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3062.                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
  3063. *
  3064. *                    Generate Q in VT
  3065. *                    (CWorkspace: need M+N, prefer M+N*NB)
  3066. *                    (RWorkspace: 0)
  3067. *
  3068.                      CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
  3069.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3070.                      IE = 1
  3071.                      ITAUQ = ITAU
  3072.                      ITAUP = ITAUQ + M
  3073.                      IWORK = ITAUP + M
  3074. *
  3075. *                    Zero out above L in A
  3076. *
  3077.                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
  3078.      $                            A( 1, 2 ), LDA )
  3079. *
  3080. *                    Bidiagonalize L in A
  3081. *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
  3082. *                    (RWorkspace: need M)
  3083. *
  3084.                      CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ),
  3085.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  3086.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3087. *
  3088. *                    Multiply right bidiagonalizing vectors in A by Q
  3089. *                    in VT
  3090. *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
  3091. *                    (RWorkspace: 0)
  3092. *
  3093.                      CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA,
  3094.      $                            WORK( ITAUP ), VT, LDVT,
  3095.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3096.                      IRWORK = IE + M
  3097. *
  3098. *                    Perform bidiagonal QR iteration, computing right
  3099. *                    singular vectors of A in VT
  3100. *                    (CWorkspace: 0)
  3101. *                    (RWorkspace: need BDSPAC)
  3102. *
  3103.                      CALL ZBDSQR( 'U', M, N, 0, 0, S, RWORK( IE ), VT,
  3104.      $                            LDVT, CDUM, 1, CDUM, 1,
  3105.      $                            RWORK( IRWORK ), INFO )
  3106. *
  3107.                   END IF
  3108. *
  3109.                ELSE IF( WNTUO ) THEN
  3110. *
  3111. *                 Path 8t(N much larger than M, JOBU='O', JOBVT='A')
  3112. *                 N right singular vectors to be computed in VT and
  3113. *                 M left singular vectors to be overwritten on A
  3114. *
  3115.                   IF( LWORK.GE.2*M*M+MAX( N+M, 3*M ) ) THEN
  3116. *
  3117. *                    Sufficient workspace for a fast algorithm
  3118. *
  3119.                      IU = 1
  3120.                      IF( LWORK.GE.WRKBL+2*LDA*M ) THEN
  3121. *
  3122. *                       WORK(IU) is LDA by M and WORK(IR) is LDA by M
  3123. *
  3124.                         LDWRKU = LDA
  3125.                         IR = IU + LDWRKU*M
  3126.                         LDWRKR = LDA
  3127.                      ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN
  3128. *
  3129. *                       WORK(IU) is LDA by M and WORK(IR) is M by M
  3130. *
  3131.                         LDWRKU = LDA
  3132.                         IR = IU + LDWRKU*M
  3133.                         LDWRKR = M
  3134.                      ELSE
  3135. *
  3136. *                       WORK(IU) is M by M and WORK(IR) is M by M
  3137. *
  3138.                         LDWRKU = M
  3139.                         IR = IU + LDWRKU*M
  3140.                         LDWRKR = M
  3141.                      END IF
  3142.                      ITAU = IR + LDWRKR*M
  3143.                      IWORK = ITAU + M
  3144. *
  3145. *                    Compute A=L*Q, copying result to VT
  3146. *                    (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
  3147. *                    (RWorkspace: 0)
  3148. *
  3149.                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
  3150.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3151.                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
  3152. *
  3153. *                    Generate Q in VT
  3154. *                    (CWorkspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB)
  3155. *                    (RWorkspace: 0)
  3156. *
  3157.                      CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
  3158.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3159. *
  3160. *                    Copy L to WORK(IU), zeroing out above it
  3161. *
  3162.                      CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ),
  3163.      $                            LDWRKU )
  3164.                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
  3165.      $                            WORK( IU+LDWRKU ), LDWRKU )
  3166.                      IE = 1
  3167.                      ITAUQ = ITAU
  3168.                      ITAUP = ITAUQ + M
  3169.                      IWORK = ITAUP + M
  3170. *
  3171. *                    Bidiagonalize L in WORK(IU), copying result to
  3172. *                    WORK(IR)
  3173. *                    (CWorkspace: need   2*M*M+3*M,
  3174. *                                 prefer 2*M*M+2*M+2*M*NB)
  3175. *                    (RWorkspace: need   M)
  3176. *
  3177.                      CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S,
  3178.      $                            RWORK( IE ), WORK( ITAUQ ),
  3179.      $                            WORK( ITAUP ), WORK( IWORK ),
  3180.      $                            LWORK-IWORK+1, IERR )
  3181.                      CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU,
  3182.      $                            WORK( IR ), LDWRKR )
  3183. *
  3184. *                    Generate right bidiagonalizing vectors in WORK(IU)
  3185. *                    (CWorkspace: need   2*M*M+3*M-1,
  3186. *                                 prefer 2*M*M+2*M+(M-1)*NB)
  3187. *                    (RWorkspace: 0)
  3188. *
  3189.                      CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
  3190.      $                            WORK( ITAUP ), WORK( IWORK ),
  3191.      $                            LWORK-IWORK+1, IERR )
  3192. *
  3193. *                    Generate left bidiagonalizing vectors in WORK(IR)
  3194. *                    (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
  3195. *                    (RWorkspace: 0)
  3196. *
  3197.                      CALL ZUNGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
  3198.      $                            WORK( ITAUQ ), WORK( IWORK ),
  3199.      $                            LWORK-IWORK+1, IERR )
  3200.                      IRWORK = IE + M
  3201. *
  3202. *                    Perform bidiagonal QR iteration, computing left
  3203. *                    singular vectors of L in WORK(IR) and computing
  3204. *                    right singular vectors of L in WORK(IU)
  3205. *                    (CWorkspace: need 2*M*M)
  3206. *                    (RWorkspace: need BDSPAC)
  3207. *
  3208.                      CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
  3209.      $                            WORK( IU ), LDWRKU, WORK( IR ),
  3210.      $                            LDWRKR, CDUM, 1, RWORK( IRWORK ),
  3211.      $                            INFO )
  3212. *
  3213. *                    Multiply right singular vectors of L in WORK(IU) by
  3214. *                    Q in VT, storing result in A
  3215. *                    (CWorkspace: need M*M)
  3216. *                    (RWorkspace: 0)
  3217. *
  3218.                      CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ),
  3219.      $                           LDWRKU, VT, LDVT, CZERO, A, LDA )
  3220. *
  3221. *                    Copy right singular vectors of A from A to VT
  3222. *
  3223.                      CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
  3224. *
  3225. *                    Copy left singular vectors of A from WORK(IR) to A
  3226. *
  3227.                      CALL ZLACPY( 'F', M, M, WORK( IR ), LDWRKR, A,
  3228.      $                            LDA )
  3229. *
  3230.                   ELSE
  3231. *
  3232. *                    Insufficient workspace for a fast algorithm
  3233. *
  3234.                      ITAU = 1
  3235.                      IWORK = ITAU + M
  3236. *
  3237. *                    Compute A=L*Q, copying result to VT
  3238. *                    (CWorkspace: need 2*M, prefer M+M*NB)
  3239. *                    (RWorkspace: 0)
  3240. *
  3241.                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
  3242.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3243.                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
  3244. *
  3245. *                    Generate Q in VT
  3246. *                    (CWorkspace: need M+N, prefer M+N*NB)
  3247. *                    (RWorkspace: 0)
  3248. *
  3249.                      CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
  3250.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3251.                      IE = 1
  3252.                      ITAUQ = ITAU
  3253.                      ITAUP = ITAUQ + M
  3254.                      IWORK = ITAUP + M
  3255. *
  3256. *                    Zero out above L in A
  3257. *
  3258.                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
  3259.      $                            A( 1, 2 ), LDA )
  3260. *
  3261. *                    Bidiagonalize L in A
  3262. *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
  3263. *                    (RWorkspace: need M)
  3264. *
  3265.                      CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ),
  3266.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  3267.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3268. *
  3269. *                    Multiply right bidiagonalizing vectors in A by Q
  3270. *                    in VT
  3271. *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
  3272. *                    (RWorkspace: 0)
  3273. *
  3274.                      CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA,
  3275.      $                            WORK( ITAUP ), VT, LDVT,
  3276.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3277. *
  3278. *                    Generate left bidiagonalizing vectors in A
  3279. *                    (CWorkspace: need 3*M, prefer 2*M+M*NB)
  3280. *                    (RWorkspace: 0)
  3281. *
  3282.                      CALL ZUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
  3283.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3284.                      IRWORK = IE + M
  3285. *
  3286. *                    Perform bidiagonal QR iteration, computing left
  3287. *                    singular vectors of A in A and computing right
  3288. *                    singular vectors of A in VT
  3289. *                    (CWorkspace: 0)
  3290. *                    (RWorkspace: need BDSPAC)
  3291. *
  3292.                      CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT,
  3293.      $                            LDVT, A, LDA, CDUM, 1,
  3294.      $                            RWORK( IRWORK ), INFO )
  3295. *
  3296.                   END IF
  3297. *
  3298.                ELSE IF( WNTUAS ) THEN
  3299. *
  3300. *                 Path 9t(N much larger than M, JOBU='S' or 'A',
  3301. *                         JOBVT='A')
  3302. *                 N right singular vectors to be computed in VT and
  3303. *                 M left singular vectors to be computed in U
  3304. *
  3305.                   IF( LWORK.GE.M*M+MAX( N+M, 3*M ) ) THEN
  3306. *
  3307. *                    Sufficient workspace for a fast algorithm
  3308. *
  3309.                      IU = 1
  3310.                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
  3311. *
  3312. *                       WORK(IU) is LDA by M
  3313. *
  3314.                         LDWRKU = LDA
  3315.                      ELSE
  3316. *
  3317. *                       WORK(IU) is M by M
  3318. *
  3319.                         LDWRKU = M
  3320.                      END IF
  3321.                      ITAU = IU + LDWRKU*M
  3322.                      IWORK = ITAU + M
  3323. *
  3324. *                    Compute A=L*Q, copying result to VT
  3325. *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
  3326. *                    (RWorkspace: 0)
  3327. *
  3328.                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
  3329.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3330.                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
  3331. *
  3332. *                    Generate Q in VT
  3333. *                    (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB)
  3334. *                    (RWorkspace: 0)
  3335. *
  3336.                      CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
  3337.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3338. *
  3339. *                    Copy L to WORK(IU), zeroing out above it
  3340. *
  3341.                      CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ),
  3342.      $                            LDWRKU )
  3343.                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
  3344.      $                            WORK( IU+LDWRKU ), LDWRKU )
  3345.                      IE = 1
  3346.                      ITAUQ = ITAU
  3347.                      ITAUP = ITAUQ + M
  3348.                      IWORK = ITAUP + M
  3349. *
  3350. *                    Bidiagonalize L in WORK(IU), copying result to U
  3351. *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
  3352. *                    (RWorkspace: need M)
  3353. *
  3354.                      CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S,
  3355.      $                            RWORK( IE ), WORK( ITAUQ ),
  3356.      $                            WORK( ITAUP ), WORK( IWORK ),
  3357.      $                            LWORK-IWORK+1, IERR )
  3358.                      CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU, U,
  3359.      $                            LDU )
  3360. *
  3361. *                    Generate right bidiagonalizing vectors in WORK(IU)
  3362. *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB)
  3363. *                    (RWorkspace: 0)
  3364. *
  3365.                      CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
  3366.      $                            WORK( ITAUP ), WORK( IWORK ),
  3367.      $                            LWORK-IWORK+1, IERR )
  3368. *
  3369. *                    Generate left bidiagonalizing vectors in U
  3370. *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
  3371. *                    (RWorkspace: 0)
  3372. *
  3373.                      CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
  3374.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3375.                      IRWORK = IE + M
  3376. *
  3377. *                    Perform bidiagonal QR iteration, computing left
  3378. *                    singular vectors of L in U and computing right
  3379. *                    singular vectors of L in WORK(IU)
  3380. *                    (CWorkspace: need M*M)
  3381. *                    (RWorkspace: need BDSPAC)
  3382. *
  3383.                      CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
  3384.      $                            WORK( IU ), LDWRKU, U, LDU, CDUM, 1,
  3385.      $                            RWORK( IRWORK ), INFO )
  3386. *
  3387. *                    Multiply right singular vectors of L in WORK(IU) by
  3388. *                    Q in VT, storing result in A
  3389. *                    (CWorkspace: need M*M)
  3390. *                    (RWorkspace: 0)
  3391. *
  3392.                      CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ),
  3393.      $                           LDWRKU, VT, LDVT, CZERO, A, LDA )
  3394. *
  3395. *                    Copy right singular vectors of A from A to VT
  3396. *
  3397.                      CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
  3398. *
  3399.                   ELSE
  3400. *
  3401. *                    Insufficient workspace for a fast algorithm
  3402. *
  3403.                      ITAU = 1
  3404.                      IWORK = ITAU + M
  3405. *
  3406. *                    Compute A=L*Q, copying result to VT
  3407. *                    (CWorkspace: need 2*M, prefer M+M*NB)
  3408. *                    (RWorkspace: 0)
  3409. *
  3410.                      CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
  3411.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3412.                      CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
  3413. *
  3414. *                    Generate Q in VT
  3415. *                    (CWorkspace: need M+N, prefer M+N*NB)
  3416. *                    (RWorkspace: 0)
  3417. *
  3418.                      CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
  3419.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3420. *
  3421. *                    Copy L to U, zeroing out above it
  3422. *
  3423.                      CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
  3424.                      CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
  3425.      $                            U( 1, 2 ), LDU )
  3426.                      IE = 1
  3427.                      ITAUQ = ITAU
  3428.                      ITAUP = ITAUQ + M
  3429.                      IWORK = ITAUP + M
  3430. *
  3431. *                    Bidiagonalize L in U
  3432. *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
  3433. *                    (RWorkspace: need M)
  3434. *
  3435.                      CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ),
  3436.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  3437.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3438. *
  3439. *                    Multiply right bidiagonalizing vectors in U by Q
  3440. *                    in VT
  3441. *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
  3442. *                    (RWorkspace: 0)
  3443. *
  3444.                      CALL ZUNMBR( 'P', 'L', 'C', M, N, M, U, LDU,
  3445.      $                            WORK( ITAUP ), VT, LDVT,
  3446.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3447. *
  3448. *                    Generate left bidiagonalizing vectors in U
  3449. *                    (CWorkspace: need 3*M, prefer 2*M+M*NB)
  3450. *                    (RWorkspace: 0)
  3451. *
  3452.                      CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
  3453.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3454.                      IRWORK = IE + M
  3455. *
  3456. *                    Perform bidiagonal QR iteration, computing left
  3457. *                    singular vectors of A in U and computing right
  3458. *                    singular vectors of A in VT
  3459. *                    (CWorkspace: 0)
  3460. *                    (RWorkspace: need BDSPAC)
  3461. *
  3462.                      CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT,
  3463.      $                            LDVT, U, LDU, CDUM, 1,
  3464.      $                            RWORK( IRWORK ), INFO )
  3465. *
  3466.                   END IF
  3467. *
  3468.                END IF
  3469. *
  3470.             END IF
  3471. *
  3472.          ELSE
  3473. *
  3474. *           N .LT. MNTHR
  3475. *
  3476. *           Path 10t(N greater than M, but not much larger)
  3477. *           Reduce to bidiagonal form without LQ decomposition
  3478. *
  3479.             IE = 1
  3480.             ITAUQ = 1
  3481.             ITAUP = ITAUQ + M
  3482.             IWORK = ITAUP + M
  3483. *
  3484. *           Bidiagonalize A
  3485. *           (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
  3486. *           (RWorkspace: M)
  3487. *
  3488.             CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
  3489.      $                   WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
  3490.      $                   IERR )
  3491.             IF( WNTUAS ) THEN
  3492. *
  3493. *              If left singular vectors desired in U, copy result to U
  3494. *              and generate left bidiagonalizing vectors in U
  3495. *              (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB)
  3496. *              (RWorkspace: 0)
  3497. *
  3498.                CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
  3499.                CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
  3500.      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
  3501.             END IF
  3502.             IF( WNTVAS ) THEN
  3503. *
  3504. *              If right singular vectors desired in VT, copy result to
  3505. *              VT and generate right bidiagonalizing vectors in VT
  3506. *              (CWorkspace: need 2*M+NRVT, prefer 2*M+NRVT*NB)
  3507. *              (RWorkspace: 0)
  3508. *
  3509.                CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
  3510.                IF( WNTVA )
  3511.      $            NRVT = N
  3512.                IF( WNTVS )
  3513.      $            NRVT = M
  3514.                CALL ZUNGBR( 'P', NRVT, N, M, VT, LDVT, WORK( ITAUP ),
  3515.      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
  3516.             END IF
  3517.             IF( WNTUO ) THEN
  3518. *
  3519. *              If left singular vectors desired in A, generate left
  3520. *              bidiagonalizing vectors in A
  3521. *              (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB)
  3522. *              (RWorkspace: 0)
  3523. *
  3524.                CALL ZUNGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ),
  3525.      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
  3526.             END IF
  3527.             IF( WNTVO ) THEN
  3528. *
  3529. *              If right singular vectors desired in A, generate right
  3530. *              bidiagonalizing vectors in A
  3531. *              (CWorkspace: need 3*M, prefer 2*M+M*NB)
  3532. *              (RWorkspace: 0)
  3533. *
  3534.                CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
  3535.      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
  3536.             END IF
  3537.             IRWORK = IE + M
  3538.             IF( WNTUAS .OR. WNTUO )
  3539.      $         NRU = M
  3540.             IF( WNTUN )
  3541.      $         NRU = 0
  3542.             IF( WNTVAS .OR. WNTVO )
  3543.      $         NCVT = N
  3544.             IF( WNTVN )
  3545.      $         NCVT = 0
  3546.             IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN
  3547. *
  3548. *              Perform bidiagonal QR iteration, if desired, computing
  3549. *              left singular vectors in U and computing right singular
  3550. *              vectors in VT
  3551. *              (CWorkspace: 0)
  3552. *              (RWorkspace: need BDSPAC)
  3553. *
  3554.                CALL ZBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), VT,
  3555.      $                      LDVT, U, LDU, CDUM, 1, RWORK( IRWORK ),
  3556.      $                      INFO )
  3557.             ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN
  3558. *
  3559. *              Perform bidiagonal QR iteration, if desired, computing
  3560. *              left singular vectors in U and computing right singular
  3561. *              vectors in A
  3562. *              (CWorkspace: 0)
  3563. *              (RWorkspace: need BDSPAC)
  3564. *
  3565.                CALL ZBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), A,
  3566.      $                      LDA, U, LDU, CDUM, 1, RWORK( IRWORK ),
  3567.      $                      INFO )
  3568.             ELSE
  3569. *
  3570. *              Perform bidiagonal QR iteration, if desired, computing
  3571. *              left singular vectors in A and computing right singular
  3572. *              vectors in VT
  3573. *              (CWorkspace: 0)
  3574. *              (RWorkspace: need BDSPAC)
  3575. *
  3576.                CALL ZBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), VT,
  3577.      $                      LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ),
  3578.      $                      INFO )
  3579.             END IF
  3580. *
  3581.          END IF
  3582. *
  3583.       END IF
  3584. *
  3585. *     Undo scaling if necessary
  3586. *
  3587.       IF( ISCL.EQ.1 ) THEN
  3588.          IF( ANRM.GT.BIGNUM )
  3589.      $      CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
  3590.      $                   IERR )
  3591.          IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM )
  3592.      $      CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1,
  3593.      $                   RWORK( IE ), MINMN, IERR )
  3594.          IF( ANRM.LT.SMLNUM )
  3595.      $      CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
  3596.      $                   IERR )
  3597.          IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM )
  3598.      $      CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1,
  3599.      $                   RWORK( IE ), MINMN, IERR )
  3600.       END IF
  3601. *
  3602. *     Return optimal workspace in WORK(1)
  3603. *
  3604.       WORK( 1 ) = MAXWRK
  3605. *
  3606.       RETURN
  3607. *
  3608. *     End of ZGESVD
  3609. *
  3610.       END
  3611.